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FOREWORD 


This report was prepared between June 1988 and September 1988 by Eagle Engineering, Inc. for 
the Advanced Programs Office of Johnson Space Center, a field center of the National Aeronau- 
tics and Space Administration. The objective is to provide a program which can be used to 
analyze the performance of a spacecraft making a low-thrust flight between the Earth and the 
Moon. 

Dr. J.W. Aired was the NASA technical monitor for the Advanced Space Transportation Study 
contract of which this task was a part. Mr. Andy Petro was the NASA task monitor for this 
particular task. Mr. W.R. Stump was the Eagle project manager. Mr. C.C. Varner was the Eagle 
task manager for this task. This program was originally written and documented in BASIC by 
Mr. D.J. Korsmeyer of the Large Scale Programs Institute at the University of Texas. The 
conversion to FORTRAN was completed by Mr. M. D’Onofrio of Eagle, and FORTRAN 
documentation was prepared by Mr. C.C Varner. 
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1.0 USING CISLUNAR 


CISLUNAR is a stand alone program designed to generate the trajectory of a low-thrust space- 
craft travelling in Earth-Moon space. The program allows the creation of functional trajectories 
dependent upon the supplied spacecraft characteristics. The trajectory generation is a user 
interactive process. The original intent was for the program user to modify the necessary control 
values until a satisfactory trajectory has been created. 

The program is started by simply typing CISLUNAR. The information that appears on the screen 
indicates that CISLUNAR has started, and shows the spacecraft’s default characteristics. These 
characteristics can be modified by the user at the beginning of each run. The program prompts 
the user for the direction of the trajectory generation by asking whether the initial orbit is about 
the Earth or the Moon. This sets the direction flags for the rest of the program. The next 
question is whether new parametric velocity curves based on the spacecraft’s characteristics 
should be created. If this question is answered "yes", the program generates these curves before 
continuing. The next three questions concern the initial altitude, velocity, and orbital position of 
the spacecraft. The altitude must be input for the program to continue. The velocity will default 
to the circular velocity at the input altitude. The final question asks if the controls for the 
spacecraft generation need to be modified. If this question is not answered or the answer is "no", 
the program uses the default values for the controls. 

Four controls are specified, Jacl, Jac2, Jac3, and Range. These four values govern the thrusting 
of the spacecraft during the final escape and translunar portions of the trajectory. Jacl indicates 
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the spacecraft is nearing the end of its spiral escape from the initial orbit; the engines shut down, 
and thrusting ceases unless the spacecraft is in the proper quadrant for transfer injection. Jac2 is 
the control that determines whether the spacecraft can achieve a cislunar trajectory. Ideally the 
Jacobian Constant at Jac2 has a value of 3±.2 <km/s>. After reaching Jac2, the spacecraft 
thrusts continuously. Jac3 is the final constraint on the amount of energy that is to be supplied to 
the spacecraft during transfer. Following Jac3, the spacecraft does not thrust. Range is the 
control that determines the distance from the initial planet that the capture guidance to the target 
planet is begun. This is the point at which reverse thrusting begins. 

Markers for these four controls show up on the trajectory as each of them is passed. For Jacl- 
Jac3, a small circle will indicate that this control has been reached. The passage of Range is 
indicated by a small vertical line. The visual representation of the controls is helpful to under- 
stand and plan a modification of the controls. The markers do not appear in FORTRAN versions 
of the program. 

While the trajectory is being generated the program can be paused, restarted, or simply stopped 
at any time. 
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2.0 PROGRAM CONSIDERATIONS 


The trajectory determination methods for impulsive and low-thrust spacecraft differ consider- 
ably. Electric propulsion systems need to thrust continuously for long periods of time in order to 
achieve a significant velocity change. Chemical, or impulsive propulsion, can create a near 
instantaneous change in the velocity. Where an impulsive thrusting spacecraft could use two 
short powerful thrusts to transfer between LEO (Low Earth Orbit) and a higher orbit, a low- 
thrust orbital transfer starting in LEO would be accomplished as a very slow outward spiral to 
the desired altitude. The fractional increase of the orbital radius per revolution is very small for 
low-thrust spacecraft. This complicates the calculation of trajectories for low-thrust vehicles. 

The characteristics of the low-thrust spacecraft will also play an important role in determining 
the type of trajectory that can be developed. The vehicle’s propulsion system and associated 
power system will have a considerable impact on the type of trajectories available. 

The power and propulsion system for a low-thrust spacecraft are intimately coupled. The 
thruster system efficiency is the fraction of electrical power that is converted to exhaust kinetic 
energy. This yields, 

m(I ip g) 2 

r\ = 

2 P 0 

{ 2 . 1 } 

where T| is the thruster system efficiency, m is the mass flow rate of the thrusters, and P 0 is the 
electrical power input to the propulsion system. 1 
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Currently, there are many different low-thrust electric propulsion systems under investigation. 
The ion engine, magnetoplasmadynamic thruster, and arcjet are a few of the leading candidates. 
All of these engines will require continuous high power to be able to perform competitively 
against chemical propulsion. Nuclear and solar power systems are the major competitors for 
high power supplies in space. Solar arrays and solar dynamic power systems have the advantage 
of using the sun as a heat source, however, they require continual sunlight. For some propulsion 
systems solar power cannot provide the needed level of power. Nuclear power, on the other 
hand, has tremendous potential for fulfilling the power needs of electric propulsion systems. 
The proposed range of power available from nuclear sources ranges from a few kilowatts to 
megawatts . 2 

This program does not include an aerocapture option. Aerocapture has numerous problems for 
large nuclear or solar power sources. The general outbound trajectory assumed is a spiral out 
from LEO and a spiral down into LLO (Low Lunar Orbit). The return trajectory is a spiral up 
from LLO and then down into LEO. 

The guidance scheme employed to determine a trajectory must use only low-thrust to capture the 
OTV into LEO. The low-thrust OTV is limited in the range of thrusting acceleration available to 
drive the vehicle to the desired orbit. Another restriction for the trajectories of nuclear-powered 
OTVs is the proposed nuclear safe orbit (NSO ). 3 This would be a designated altitude below 
which the nuclear powered spacecraft would be prohibited. The spacecraft would be prohibited 
from descending below this altitude at any point of the trajectory. 
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In the development of trajectories for low-thrust cislunar OTVs, little attention has been directed 
at the guidance and control of the spacecraft. The premise that the guidance of the vehicle and 
the determination of the appropriate trajectory are unrelated is false. Rather guidance and 
trajectory determination for low thrust vehicles are closely related problems which, by necessity, 
must be treated with equal importance. 4 


1. Hill, P.G., and Peterson, C.R., p. 336. 

2. English, Robert E., "Power Generation from Nuclear Reactors in Aerospace Applica- 
tions, "NRC Symposium on Advanced Compact Reactors, Washington, D.C., November 
15-17, 1982. 

3. Galecki, Diane L., and Patterson, Michael J., "Nuclear Powered Mars Cargo Transport 
Mission Utilizing Ion Propulsion," AIAA/SAE/ASME/ASEE 23rd Joint Propulsion 
Conference, San Diego, CA, 1987, AIAA-87-1903. 

4. Battin, R.H., and Miller, J.S., "Trajectories and Guidance Theory for a Continuous Low- 
Thrust Lunar Reconnaissance Vehicle," 6th Symposium on Ballistic Missile and 
Aerospace Technology, 1961. 
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3.0 PROBLEM FORMULATION 


A major problem in the design of low-thrust OTVs and their associated trajectories is the lack of 
an end to end simulation tool for the spacecraft trajectory, from NSO travelling to LLO and the 
subsequent return. The current concern is how the vehicle will behave at the proposed thrust 
level and how it will be guided on its trajectory. 

To adequately understand the dynamics of motion of the low-thrust spacecraft, the gravitational 
effects of the Earth and the Moon on the spacecraft must be included for the full duration of the 
trajectory. The thrusting acceleration for low-thrust OTVs in high Earth orbit is the same 
magnitude as the perturbing force due to the Moon. To model the Earth-Moon system with the 
necessary accuracy and achieve computational efficiency, the restricted three-body formulation 
of the dynamical equations is utilized as the governing equations of motion. 

3.1 RESTRICTED THREE-BODY FORMULATION 

The problem of three bodies was first formulated in 1772 by Lagrange. Further studies by 
Poincare, Laplace, Hill and Szebehely have resulted in a detailed treatment of the problem and a 
general understanding of the interactions between the two primary gravitational fields. Various 
formulations are available to represent the problem of three bodies. The formulation used in this 
study was referenced from Kaplan 1 and Moulton. 2 
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Many realistic orbital cases may be modelled as restricted three-body situations. An exemplary 
case is that of a spacecraft moving in the Earth-Moon system. Certain assumptions are made 
about the nature of the Earth-Moon system that permit a straightforward solution at a slight loss 
of accuracy. The motion of the Earth-Moon system is assumed to be circular and coplanar about 
its center of mass (barycenter) and the spacecraft, at point P, has negligible mass. This system is 
shown in Figure 3.1. The motion of the spacecraft is governed by the relative gravitational 
attraction of the Earth and the Moon rotating about the barycenter. The spacecraft is assumed to 
have no impact on the motion of the Earth or the Moon. Thus, the acceleration at P is 





{3.1} 


The absolute acceleration of the spacecraft is obtained in terms of the rotating coordinate system, 
x,y,z, by relating the acceleration of the spacecraft in the non-inertial (barycenter) rotating 
system to that in the inertial coordinate system. Hence, 


a p = a c + n x (n x r) + r b + 2n x r b 


where, 

r is the radius vector of the spacecraft, 

r k is the apparent velocity of the spacecraft in the rotating coordinates, 
r b is the apparent acceleration of the spacecraft in the rotating coordinates, 
a p is the acceleration of the spacecraft in inertial coordinates. 


{3.2} 
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(cjn.) 


jdy Problem Nomenclature 


a„ is the acceleration of the origin in the inertial coordinates, 
n is the angular velocity vector of the Earth-Moon system, n = ni^, 
n x (n x r) is the centrifugal acceleration, and 


2 n x r b is the Coriolis acceleration due to the motion of the spacecraft in x,y,z. 


Equating {3.1 } and (3.2), noting that the acceleration of the origin is the same for both equa- 
tions, and expressing the acceleration in component form, yields the equations of motion for the 
spacecraft in the rotating coordinate system. 


x - 2ny - n 2 x = 


y + 2nx - n 2 y = 



{3.4} 


3.2 JACOBIAN CONSTANT 


In this formulation of the equations of motion, the energy of the spacecraft is not conserved. 
However, the sum of the angular momentum, velocity, and potential energy of the spacecraft is 
conserved. This can be shown with the Jacobian Integral. Multiplying the first equation of 
(3.4) by dx/dt, the second by dy/dt, adding, and integrating the result yields this integral. 


x 2 + y 2 - n 2 (x 2 + y 2 ) = 


2 ^ 

r. + r„ - C 


{3.5} 


where C is known as the Jacobian constant. Mathematician Karl Gustav Jacobi first formulated 
this integral in 1836. This constant, C, can be determined for any set of initial conditions. 
Equation 3.5 determines the locus of those points where the spacecraft can travel given the initial 
conditions. In particular, if the velocity of the vehicle is set equal to zero for a given C, equation 
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3.5 will describe a curve where the spacecraft’s motion is bounded. On this curve the spacecraft 
with a given C will have zero velocity. Only on the inside of the curve will the square of the 
spacecraft’s velocity be positive, restricting the motion of the vehicle to that side. Figure 3.2 
shows a series of zero velocity curves in the Earth-Moon system. 


1. Kaplan, Marshall H., Page 290-300. 

2. Moulton, Forest Ray, An Introduction to Celestial Mechanics . (Dover Publications, Inc., 
New York: 1914, 2nd Revised Edition) pages 277-287. 
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Figure 3.2 Zero Velocity Curves in the Earth-Moon System 
(Kaplan, Page 292) 


11 


4.0 GUIDANCE METHODOLOGY 


The trajectory determination and guidance of a cislunar low-thrust OTV is divided into three 
distant phases: departure, translunar targeting, and capture. Each of these phases has a different 
guidance scheme to achieve the overall goal of generating a trajectory between the Earth and the 
Moon. 


4.1 DEPARTURE 

The first phase in any cislunar journey for an OTV is the escape from the initial parking orbit, 
whether about the Earth or the Moon. For low-thrust spacecraft to achieve escape, a long period 
of continuous thrusting is necessary. This results in a slowly increasing spiral trajectory from 
the initial orbit. The direction of the thrust vector should be in the direction that has the highest 
rate of increase of the energy of the orbit per revolution. It can be shown that a near-optimal 
thrust for an orbital transfer should be directed along the velocity vector of the spacecraft for the 
majority of the trajectory. 12 This is referred to as tangential thrust, because the thrusting 
acceleration will be tangent to the trajectory at all times. Another thrusting scheme, circum- 
ferential thrust, directs the acceleration along a vector perpendicular to the central body. 
Tangential thrust is the thrusting approach used in the spiral escape from the departure planet. 


12 



4.2 TRANSLUNAR TARGETING 


The value of the Jacobian constant of a spacecraft will be used as an indicator of the sufficient 
energy for the cislunar transfer. Dr. Victor Szebehely 3 notes that an equipotential curve like that 
shown in Figure 4.1 occurs when the Jacobian constant is approximately 3.3. When the space- 
craft is outside of this curve, the range of motion is no longer restricted only to geocentric or 
lunar orbit, but can include transfer between the neighboihood of the Earth or Moon. 

The spacecraft needs to achieve the required Jacobian constant when the velocity vector of the 
spacecraft is pointed in the appropriate direction to allow transfer between the Earth and the 
Moon. Figure 4.2 shows this targeting procedure for the OTV. The area about the Earth is 
divided into four quadrants, I-IV. The Jacobian of the spacecraft is calculated continuously as 
the spacecraft nears escape. Various values of the Jacobian are chosen experimentally to act as 
indicators of the spacecraft’s proximity to escape. The initial indicator of escape is the value of 
the Jacobian while the spacecraft is in the third quadrant. If the vehicle achieves a Jacobian of 
4.1 while located in the third quadrant, continued thrust will enable a lunar passage to occur. 
However, if the spacecraft achieves the value of the initial Jacobian, 4.1, while outside the third 
quadrant, the spacecraft’s thrust is turned off. When the spacecraft arrives in the third quadrant 
the thrust is reinitiated tangentially to obtain the necessary Jacobian for escape and remains on 
until the spacecraft achieves sufficient energy for transfer and enters the capture phase. On the 
return trip back from the Moon, the same methodology is used with different Jacobian constants. 
The Jacobian constants used in the Earth to Moon voyage are driven only by the acceleration 
level of the spacecraft during escape. 
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4.3 CAPTURE 


As the vehicle approaches the Moon the capture guidance phase of the trajectory is initiated. In 
the absence of impulsive thrust, the approach and capture to the target body are critical and must 
not necessitate maneuvering beyond the limited capabilities of the propulsion system. The 
problem of low-thrust spacecraft guidance and trajectory determination between the Earth and 
the Moon was addressed in a study by Richard H. Battin and James S. Miller in the late 1950’s 
and early 1960’s. The concept for the spacecraft guidance during capture used in this study is 
derived from Battin and Miller’s work. 4 

The guidance scheme is relatively simple and straightforward. The operation of the capture 
phase guidance is illustrated in Figure 4.3. The velocity of the spacecraft, V y , relative to the 
target body (i.e. the Earth or the Moon) is compared with a precalculated velocity, V e , profile for 
a spiral capture. This velocity profile is a function of the radial distance from the target body 
and the magnitude of the thrust acceleration. This velocity difference, V d , is used in combination 
with the nominal acceleration to determine the direction and magnitude of the spacecraft thrust 
during capture. 

In order to calculate the velocity as a function of the radial distance from the capturing body, the 
"ideal", or reference trajectory must be calculated. This is a spiral capture that achieves circular 
velocity at the desired final altitude. To determine this reference spiral path and the velocity 
vectors that accompany it, the spacecraft starts in a circular orbit at the desired final altitude 
about the target body. The mass of the spacecraft at the final altitude is determined by 
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Jacobian of 4.1 



Figure 4.2 - Translunar Targeting 
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Figure 4.3 - Capture Phase Thrusting Guidance 
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estimating the final mass of the spacecraft at the completion of the mission to be 80% of the 
initial mass. The spacecraft follows a spiral out from the target planet using tangentially directed 
thrust and a negative mass flow. Only the gravitational field of the target planet is considered. 
The spiral trajectory is otherwise without perturbations and consequently remains two-dimen- 
sional. The calculation of the trajectory continues until the energy of the orbit is non-negative, 
and the vehicle is on a parabolic path. The associated radial and tangential components of the 
spacecraft’s velocity are recorded at select steps as functions of the radial distance from the 
central body. The velocity functions are obtained by fitting the recorded velocity components to 
polynomial and power curves. Figures 4.4 and 4.5 are example graphs of the velocity profiles 
for tangential and radial velocity at a radial distance. The equations shown in the figures have 
parameterized the velocities as a function of the radial distance. This data was obtained by the 
described reverse integration process. 

An explanation of the thrust guidance control used by the spacecraft during capture phase of the 
trajectory is presented as follows. The actual velocity of the spacecraft, V„ at a given radial 
distance, r, is compared with the parameterized reference capture velocity, V c , at r. The 
difference between these velocity vectors is then determined as V d , where 


The instantaneous change in the capture velocity can be approximated as the effect of the 
acceleration of the spacecraft due to its thrust and the gravitational pull of the planet acting over 


a small time increment, At. This implies 
V c » <a c + g) At 
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Tang. Velocity (km/s) 


Tangential Velocity .vs. Radial Distance 



Figure 4.4 Tangential Velocity as a Function of Radial Distance 
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Radial Velocity (km/s) 


Radial Velocity .vs. Radial Distance 



Figure 4.5 Radial Velocity as a Function of Radial Distance 
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where a,, is the nominal acceleration of the spacecraft, and g is the gravitational acceleration 
vector of the capture planet. The instantaneous change in the actual velocity of the spacecraft 
can also be approximated as such 

V v = (a t + g) At 


{4.3} 


where a, is the acceleration vector of the spacecraft on the trajectory. Combining equations 4.1, 
4.2, and 4.3 and rearranging to find a t , equation 4.4 is obtained. 


Av d 



The thrust acceleration is then chosen so that the rate of change of the velocity vector V d is 
proportional to V d itself. This results in 
Av d AV d 

At T c 

where T c is an empirically determined time constant. With this formula the appropriate thrust 
acceleration can be determined in both magnitude and direction simply with the knowledge of 
the vehicle’s position, velocity, and nominal thrust acceleration a c . 


In the application of the guidance it is reasonable to assume that the direction of the thrust 
acceleration can be varied at will, but the magnitude of the thrust is limited by the capabilities of 
the propulsion system. The spacecraft thrust is never required to deliver greater than the 
nominal thrust. The possibility of a reduction in the thrusting acceleration is not precluded as a 
desirable effect of the thrusting algorithm. Figure 4.6 is a graphical representation of the 
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acceleration vectors a, and a c . The radii of the circles are determined by the nominal acceleration 
of the spacecraft. 

Then, 

z* 

«t < a= + T c 

When the magnitude of the thrust acceleration, a,, is less than the nominal capabilities of the 
engine, a less than nominal thrust is needed in the appropriate direction. 

This thrusting algorithm drives the spacecraft’s velocity components toward the reference 
velocity conditions. The actual reference conditions can never be reached due to the fact that 
they were generated in an ideal two-body environment with estimated final conditions, and most 
importantly, the magnitude of the spacecraft’s thrust is limited. This results in a generated spiral 
capture trajectory that is not "ideal" but is adequate. 


1. Keaton, Paul W., page 4 

2. Hill, P.G., and Peterson, C.R., Chapter 10. 

3. Szebehely, Victor, Personal Communication, November, 1987. 

4. Battin, R.H., and Miller, J.S. 


22 


i 




23 


5.0 PROGRAM INPUTS 


Cislunar needs some preliminary data before it can operate. This data is provided by the user 
interactively during program execution. A series of screen prompts will direct the user to supply 
vital information. The user will input the requested data; and the program will proceed to the 
next question. In this section the screen prompts are discussed as well as the information that the 
program expects the user to supply. 


1 . Prompt: DO YOU WISH TO SUPPLY S/C CHARACTERISTICS ( Y OR N) 

Description: Type "Y" if the default values are to be changed. After hitting RETURN the 
program will request further information about the spacecraft characteristics. 
Typing "N" at this prompt instructs the program to use default spacecraft 
characteristics. 


1 A. Prompt: 
Description: 


SPACECRAFT INITIAL MASS = 

Enter the spacecraft mass before it leaves the holding orbit about the planet of 
origin. The mass is in kilograms. The default value is 60,000 kg. 


IB. Prompt: 
Description: 


SPECIFIC IMPULSE OF ENGINE = 

Enter the engine specific impulse in seconds. Default value is 2500s. 
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1C. Prompt: 


MASS FLOW RATE OF ENGINES = 


Description: 

ID. Prompt: 
Description: 


2. Prompt: 
Description: 

3. Prompt: 
Description: 


Input the propellant mass flow rate. The units are in kilograms per second, 
and the default value is 0.0082 kg/s. 

DEGREE OF POLYNOMIAL CURVE FIT (2-7) 

When requested the program creates a velocity guidance spiral about the 
target planet. This guidance spiral is used by the control system to provide 
capture targeting. After creating the data for this guidance spiral, the program 
forms a curve which approximately "fits" the data. This curve is described 
with a polynomial mathematical expression. The expression can be between 
2nd through 7th order; and the order must be supplied by the user. Note: 
Higher order curves take longer to create than lower order curves. 3rd order 
is usually sufficient. 

STARTING ORBIT ABOUT THE EARTH OR THE MOON? (E OR M) 
Type the first letter of the planetary body from which the spacecraft origi- 
nates. 


DO YOU WANT TO GENERATE PARAMETRIC VELOCITY CURVES? 
If the user types "Y" or "Yes" then the program creates the velocity guidance 
spiral discussed in question ID. In FORTRAN versions of the program, the 
parametric velocity curves are automatically generated; this question is not 
asked. In the BASIC version of CISLUNAR, the user is given the option due 
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4. Prompt: 
Description: 

5. Prompt: 
Description: 


6. Prompt: 
Description: 


to the length of time required to generate these curves. The initial attempts at 
cislunar targeting are not likely to come close enough to the target planet to 
make capture guidance worthwhile. In such cases, the velocity guidance 
spiral generation is not only unnecessary, but also time consuming. If the 
parametric velocity curves are not desired, the user should type "N" or "No" 
at this prompt. 

INPUT THE RADIUS FROM THE PLANET’S CENTER 

Enter the radial distance of the spacecraft from the center of the planet of 

origin. This distance has units of kilometers. 

INPUT ANGLE (DEG.) FROM - X AXIS 

The X axis is the line between the Earth and the Moon. The -X axis is the 
Earth-Moon line on the "Moon" side of the Earth, or on the "Far" side of the 
Moon. While the +X axis is the Earth-Moon line on the "Far" side of the 
Earth, or the "Earth" side of the Moon. The program requests that the user 
supply the angle, measured counter-clockwise, from the -X axis. The angle is 
in units of degrees. 

DO YOU WISH TO SPECIFY THE VELOCITY? (Y OR N) 

Enter "N" to tell the program to default to circular orbit speed, otherwise enter 
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6A. Prompt: INPUT ABSOLUTE VELOCITY 

Description: Enter spacecraft velocity in kilometers per second. 

7. Prompt: MODIFY CONTROL JACOBIANS AND RANGE? 

Description: The Jacobian Constant is the controlling parameter for cislunar targeting. The 
RANGE parameter is used to initiate capture guidance and thrust control. The 
user should type "N" at this prompt if the default values are suitable. 
Otherwise, enter "Y" and proceed to modify the control Jacobians and the 

RANGE. 

1 A. Prompt: JAC 1 = 4.93 

Description: Input the new value for the first Jacobian control point. The Jacobian 
constant is a non-dimensional number, which can be represented by zero- 
velocity carves as shown in Figure 4.1. After passing the first control point 
the spacecraft will thrust only while in control quadrant. 

7B. Prompt: JAC 2 = 4.10 

Description: Enter the new value for the second Jacobian control point. The passage of the 
second control point means that the thrusting will revert back to continuous 
mode. 
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7C. Prompt: 


JAC 3 = 2.70 


Description: 


7D. Prompt: 
Description: 


Enter the new value for the third Jacobian control point. The engines shut 
down entirely after reaching the third control point. Thrust is zero and the 
spacecraft coasts. 


RANGE = 255000 

Input the new RANGE in kilometers. The RANGE is the distance from the 
planet of origin at which capture guidance is to begin. 
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6.0 TEST CASE 


An example of a set of inputs for cislunar flight from the Eaith which works well is the 
following: 


Prompt: DO YOU WISH TO SUPPLY S/C CHARACTERISTICS? (Y OR N) 

Answer: "N" 

Prompt: STARTING ORBIT ABOUT THE EARTH OR MOON? (E OR M) 
Answer: "E" 

Prompt: DO YOU WANT TO GENERATE PARAMETRIC VELOCITY CURVES? 
Answer: "Y" 

Prompt: INPUT THE RADIUS FROM THE PLANET’S CENTER 
Answer: 19000 

Prompt: INPUT ANGLE (DEG) FROM -X AXIS 
Answer: -76 
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Prompt: DO YOU WISH TO SPECIFY THE VELOCITY? (Y OR N) 

Answer: "N" 

Prompt: MODIFY CONTROL JACOBIANS AND RANGE? 

Answer: "Y" 

Prompt: JAC 1 = 4.93 
Answer: 4.93 

Prompt: JAC 2 = 4. 10 
Answer: 4.10 

Prompt: JAC 3 = 2.70 
Answer: 2.55 

Prompt: RANGE = 255000 
Answer: 255000 

The graphical results are shown in Figures 6.1 and 6.2. The run data is stored in a file called 
CISLUNAR.OUT. This file is not created in BASIC versions since the data is to be presented 
with the graphics directly on the screen. 

The total velocity changes that a vehicle must undergo to perform this mission are derived from 
Tsiolkovsky’s equation (6.1). 
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Figure 6.1 - CISLUNAR Output: BASIC Version 



Figure 6.2 - CIS LUNAR Output: FORTRAN Version 
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AV = g * I, p * LN 


{ 6 . 1 } 


[ ] 

M “ 

Av = Total Velocity Change <km/s> 
g = Surface gravity of the Earth (9.81 x 10~ 3 km/s) 
I jp = Spacecraft's Specific Impulse <s> 

M 0 = Spacecraft's Initial Mass <kg> 

M = Spacecraft's Final Mass <kg> 

AV = 4.619 <km/s> 


In order to check the validity of this program, let us compare these results with those that we 
would receive from other methods of delta-V determination. One method of approximating the 
total delta-V required for very-low thrust orbital transfers is described in the following para- 
graph, and will hereafter be called the Approximation of Low-Thrust Velocity Change. 


The hypothesis behind the Approximation of Low-Thrust Velocity Change is that the low-thrust 
delta-V required to transfer between two orbits of a central force body is approximately equal to 
the difference in their mean orbital speeds. For circular orbits, this is the difference between the 
circular orbit speeds of the two orbits between which the vehicle is to transfer (Equation 6.2). 


Av LowThrujt = Absolute Value of (V^ - V^) 

{ 6 . 2 } 

where: V ml = Mean Orbital Velocity of Orbit #1 

V^ = Mean Orbital Velocity of Orbit #2 


In the test case shown, the spacecraft makes a low-thrust transfer about the Earth between the 
orbits of 19,000 km circular and the Moon’s Orbit of 384,400 km circular. It then makes a low- 
thrust transfer about the Moon from Escape orbit (V m = 0) to 5,000 km circular. The low-thrust 
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delta-V for these transfers are calculated using Equation 6.2 to be 3.570 <km/s> and 0.990 
<km/s> for the Earth and Moon transfers respectively. Therefore, using the Approximation of 
Low-Thrust Velocity Change, the total low-thrust delta-V is approximately 4.560 <km/s>. This 
compares well with the 4.619 <km/s> delta-V, obtained using the mass ratio, and Tsiolkovsky’s 
equation. 
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7.0 PROGRAM CODING 


The discussion of the program coding is sectioned by subroutine. The subroutines are addressed 
in alphabetical order following the discussion of the main "driver" routine. Each subroutine has 
a description, a list of variables, followed by the actual program code. There are two sets of 
code. The first set is the code for the BASIC Version of the CISLUNAR; the second set is the 
FORTRAN Version. 
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Subroutine Description 
And Variable Dictionary 


MAIN PROGRAM DESCRIPTION 


The main program controls the flow of the entire program. Initially it sets the global constants 
and the global variables. It calls the 10 subroutine which then returns the information necessary 
to begin the trajectory generation. An integration loop is run to generate the trajectory. While in 
the loop, the program determines where along the trajectory the spacecraft is and adjusts the 
integration step size for accuracy and convenience. During the translunar portion of the trajecto- 
ry the instantaneous Jacobian constant is calculated based on the spacecraft’s position and 
velocity by calling JACOBI. The main program compares this Jacobian value with the chosen 
control variables and marks an indicator along the trajectory when the controls are reached. The 
spacecraft’s mass is decremented according to the integration step size and the mass flow rate, 
and a new acceleration level for the next pass through the integration loop is calculated. The 
fourth order Runge-Kutta routine is called to integrate the position and velocity of the spacecraft. 
The position, velocity, mass, and elapsed time of flight output is updated every fifth integration. 
In the BASIC Version, each update is sent to the screen. In the FORTRAN Version, each update 
is stored in two arrays called GRAFX and GRAFY. These arrays are plotted at the end of the 
simulation. At the end of each integration loop the main program checks to see if there has been 
input from the keyboard to pause, continue, restart, or quit. 
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PROGRAM WIDE COMMON VARIABLES 


ACCEL 1 


Isp 

Mdot 

MU 

RANGE 

SCMASS 

SCMASSV 

TH 

THRUST 

VRN 

VTN 


Acceleration of spacecraft during the spiral reference velocity parametrization 
portion (km/s A 2). 

Specific Impulse of the propulsion system (seconds). 

Mass flow rate of the propulsion system (kg/s). 

Gravitational parameter of the destination planet (km A 3/s A 2). 

Range from departure planet that the capture phase is initiated (km). 

Initial spacecraft mass (kg). 

Instantaneous spacecraft mass, accounting for the propellant used (kg). 

Numeric indicator of the spacecraft’s thrust, on (I) or off (0). 

Thrust of the propulsion system (kg*km/s A 2). 

Degree of the radial velocity polynomial curve fit. 

Qass of equation for tangential velocity parameterization, linear (1), exponential 
(2), power (3), or logarithmic (4). 


PROGRAM WIDE CONSTANTS 


DE Distance of the Earth’s center from the Barycenter of the Earth-Moon system 

(km). 

DM Distance of the Moon’s center from the Barycenter of the Earth-Moon system 

(km). 

EMDIST Distance between the center of the Earth and the center of the Moon (km), 
gravity Earth’s surface gravity (kin/s A 2). 

MUE Gravitational parameter of the Earth (km A 3/s A 2). 

MUM Gravitational parameter of the Moon (km A 3/s A 2). 

MUN Ratio of the Moon’s mass to the mass of the Earth-Moon system 
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NUM 


Order of the X state vector for the Runge-Kutta. 

NW Mean motion of the Earth-Moon system (rad/s). 

NW2 Mean motion squared (rad/s) A 2. 

PI 71 

ADDITIONAL MAIN PROGRAM VARIABLES 

ACCEL Acceleration of the spacecraft during the cislunar trajectory (km/s A 2). 

ALPHA Theta + PI/2, angle of the tangential velocity vector. 

CJ Instantaneous Jacobian constant of the spacecraft. 

counter Counter for determining when to update the screen display of the trajectory. 

DIRECTION Numeric indicator of the direction of the trajectory generation. Earth to Moon (0) 
or Moon to Earth (1). 


dt Integration step size (seconds). 

RANGEOFF Flag for GUIDE subprogram indicating if the spacecraft has entered the capture 




phase of the trajectory generation. 

— 

RJAC() 

The array of the Jacobian constants used as controls for the departure portion of 
the trajectory generation. 

— 

RMT 

Radial distance of the spacecraft from the center of the Moon (km). 


RR 

Radial distance of the spacecraft from the center of the Earth (km). 


test$ 

String variable for controlling program from the keyboard. 


theta 

Angle between the radial vector to the spacecraft from the controlling gravitation- 
al body and the x-axis. Dependent upon xf. 

— 

thrst 

Flag for the spacecraft during the spiral escape indicating the passage of the 
second control Jacobian, (0) off (1) on. 

— 

TIC 

Flag for graphically showing the position on the trajectory where the control 
Jacobians and Range are reached. 


TT 

Total time of trajectory generation (seconds). 
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VR() 

VT() 

VY 

X() 

xf 


The array of parameterized radial velocity component coefficients. 

The array of the parameterized tangential velocity component coefficients. 
Magnitude of the spacecraft’s velocity (km/s). 

State Vector of the spacecraft’s position and velocity in the rotating x,y coordi- 
nates (km and km/s). 

Distance along the jt-axis the spacecraft is from the Earth or Moon. Dependent 
upon the phase of the trajectory generation, departure or capture. 
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Cislunar Program Flowchart 
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BASIC CODE 
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DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 

DECLARE 


SUB 10 (VR#(), VT#(), X#(), DIRECTION#, RJAC#<)) 

SUB DERI (X# ( ) , DX# () ) ’ 

SUB RUK (X# () , N!, dt#) 

SUB SPIRAL <X#<), DIRECTION#, RANGE SC#, VT# () , VR# ()) 
FUNCTION ACOS# (X#) u# 

SUB POLYFIT (PR# () , PV#(), DEGREE#, N#, VR# () ) 

SUB CUKVJ (X# () , Y#(), N #, VT# () ) 

SUB DBOX (urow%, ucol%, lrow%, lcol%) 

FUNCTION ATAN2# (X#, Y#) 

SUB SHOW (X#()) 

SUB JACOBI (X# () , CJ#) 

SUB RUK4 (X# () , N!, dt#) 

SUB DER (X# () , DX# () ) 

SUB GUIDE (GDI#, GD2#, X#()) 


Main: 


^Main^Prograin f° r Low-thrust Guidance; 3-body, eqns from Kaplan 


DEFDBL A-H, K-Z 

COMMON SHARED TH, SCMASSV, SCMASS, 
COMMON SHARED ACCEL1, MU, VTN, VRN 
REM Constants 
CONST MUE - 398600.5# 

CONST MUM - 4902.794# 

CONST DE - 4670.6778# 

CONST DM - 379729.32# 

CONST EMDIST - DE + DM 
CONST NW - .000002665314572# 

CONST NW2 - NW * NW 

CONST gravity - 9 . 809999999999999D-03 

CONST PI - 3.141592654# 

CONST NUM - 4 ! 

CONST MUN - MUM / (MUM + MUE) 

SCMASS - 60000# 

Isp - 2500# 

Mdot - 8 . 2000 00000000 00 1D-0 3 

VRN - 7 

DIM X ( 4 ) , VR ( 0 TO 7), VT(2), RJAC(3) 
AGAIN: 

IO VR(), VT ( ) , X ( ) , DIRECTION, RJAC() 
SCMASSV - SCMASS 

TT - 0#: counter - 5: thrst •» 0: TIC 

DO 


Isp, THRUST, Mdot, RANGE 


'mu of the earth 
'mu of the moon 

'distance of Earth from barycenter (Jan) 
'distance of Moon from barycenter (km) 
'distance between the earth and the moon 
'mean motion (rad/s) 

'mean motion squared 

'earth's surface gravity (km/s''2) 

'Pi 

' order of state vector for Runge-Kutta 
' ra tio of Moon's mass to system mass 
'total s/c mass (kg) 

'specific impulse (secs) 

'mass flow rate (kg/s) 

'degree of polynomial curve fit 

'program position at a restart 


1: RANGEOFF - 0' initial time 


W - SQR (X (3) A 2 + X (4 ) - 2) 

RR “ SQR ( (X (1 ) + DE) A 2 + X (2 ) ~ 2) 

RMT - SQR ( (X ( 1 ) - DM) ~ 2 + X (2) ~ 2) 

IF X(l) < 210000 THEN 'earth portion 

xf - X(l) + DE 

dt - INT(.01# * RR) 'integration step 
IF RR < 7800 AND DIRECTION - 1 THEN dt - 
IF RR > 40000 AND DIRECTION - 0 THEN 
JACOBI X(), CJ 
LOCATE 5, 13: PRINT CJ 

ELSE 


size 

5# 


CJ - 10 

END IF 

ELSE 

xf - X(l) - DM 
dt - I NT ( . 02# * RMT) 

IF RMT < 2200 AND DIRECTION 0 THEN dt - 3# 
IF RMT > 6000 AND DIRECTION ■ 1 THEN 
JACOBI X(), CJ 
LOCATE 5, 13: PRINT CJ 

ELSE 


'lunar portion 
' integration steD si?P 


END IF 


CJ 

END IF 
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ORIGINAL PAGE IS 
OE POOR QUALITY 


TIC - 2 
TIC - 3 
TIC - 4 


" kg" 


SELECT CASE TIC 

CASE 1: IF CJ < RJAC(TIC) THEN CIRCLE <X<1), 

CASE 2: IF CJ < RJAC(TIC) THEN CIRCLE (X(l), 

CASE 3: IF CJ < RJAC(TIC) THEN CIRCLE (X<1), 

CASE ELSE: TIC - 4 
END SELECT 


X(2)), 2000: 
X<2)), 2000: 
X (2) ) , 2000: 


theta - ATAN2 (xf , X(2)) 

ALPHA - theta + PI / 2# 

TT - TT + dt 

SCMASSV - SCMASSV - Mdot * dt 
ACCEL “ THRUST / SCMASSV 
RUK4 X() , NUM, dt 
IF counter >«■ 5 THEN 


'angle of radius vector 
'angle of tangential vector 

* TH' True s/c mass as propellant is used 
'acceleration of the s/c (km/s~2) 


LOCATE 3, 
LOCATE 3, 
LOCATE 3, 
LOCATE 2, 
LOCATE 2, 


17: 

35: 

67: 

67: 

40: 


PRINT 

PRINT 

PRINT 

PRINT 

PRINT 


USING 

USING 

USING 

USING 

USING 


"♦###.♦♦"; TT / 3600#; : PRINT ■ 
"#.####"; W; : PRINT " kro/s" 
•#♦#♦##.##"; RMT; : PRINT " km" 
"######.##"; RR; : PRINT " km" 
"#***#.###"; SCMASS - SCMASSV; : 


hrs~ 


pri: 


counter ~ 0 
SHOW X() 

ELSE 

counter - counter + 1 

END IF 

testS - UCASES (INKEYS) 

IF testS “ "P" THEN 
DO 


TestlS - UCASES (INKEYS) 
LOOP UNTIL TestlS - "G" 

END IF 

IF testS - "R" THEN GOTO AGAIN 
LOOP UNTIL testS - "Q" 

CLS 

END 
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FORTRAN CODE 
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1 C====== - ^== = = ========== 

2 c Main Program for Low-thrust Guidance; 3-body, eqns from Kaplan 

3 C— ======= : ====== :== — ====— ====== — “ ===== — — ==== — :==== 

4 C PROGRAM WRITTEN BY DAVID KORSMEYER 

5 C 

6 C NASA CONTRACT NAS 17878 

7 C 

8 C ELECTRIC PROPULSION; TO 87-57; TASK 1 .3 

9 C PROGRAM TRANSLATED AND MODIFIED BY 

10 C MIKE D’ONOFRIO AND CHRIS VARNER 

1 1 C EAGLE ENGINEERING, INC. 

12 C AUGUST 10,1988 

13 C 

14 IMPLICIT REAL* 16(A-H,K-Z) 

15 REAL* 16 Isp, MALT 

16 REAL*4 XPLOT,YPLOT,GRAFX,GRAFY 

17 INTEGER VRN, counter, DIRECT, NUMAU,TH,thrst,VTN 

18 CHARACTER*24 TITLED 

19 DIMENSION X(4), VR(10), VT(5), RJAC(3), GRAFX(20000), 

20 * GRAFY(20000), XPLOT(201), YPLOT(201) 

21 C 

22 C Open Graphics Routines 

23 C 

24 CALL JBEGIN 

25 CALL JDINIT ( 1 ) 

26 CALL JDEVON ( 1 ) 

27 CALL JIENAB ( 1, 4, 1 ) 

28 C 

29 C OPEN FILE FOR OUTPUT 

30 C 

3 1 OPEN (UNIT = 1 , FILE = ’CISLUNAR.OUT’ , STATUS = ’NEW’) 

32 C 

33 C INITIALIZE VARIABLES 

34 C 

35 TH = 0 

36 MUE = 398600.5 

37 MUM = 4902.794 

38 DE = 4670.6778 

39 DM = 379729.32 

40 EMDIST = DE + DM 

41 NW=.0000026653 14572 

42 NW2 = NW * NW 

43 gravity = 9.809999999999999D-03 

44 PI= 3.1415926535 

45 NUM = 4 

46 MUN = MUM / (MUM + MUE) 

47 SCMASS = 60000.0 
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48 TMAX = 600.0 * 3600.0 

49 Isp = 2500.0 

50 Mdot = 8 .20000000000000 1 D-03 

51 ACCEL1 = 0.0 

52 VRN = 7 

53 CALL IO (S CM ASS, Isp, Mdot, VRN, THRUST, X, gravity, 

54 + TITLED, RANGE, DIRECI\MUE, VR, rl, theta, RJAC, 

55 + PI, VT,VRR,MU,ACCEL1 ,MUM,NUM,DM,DE,dt,STPALT,TT,VTN,VRO) 

56 SCMASSV = SCMASS 

57 TT = 0.0 

58 counter = 5 

59 thrst = 0 

60 RANGEOFF = 0 

61 RR = QSQRT((X(1) + DE)**2. + X(2)**2. ) 

62 RMT = QSQRT((X(1) - DM)**2. + X(2)**2. ) 

63 PRINT *, ’ PLEASE WAIT, GENERATING GRAPH AND DATA FILE’ 

64 DO WHILE( DIRECT .EQ. 0 AND. RMT .GT. STPALT 

65 * .OR. DIRECT .EQ. 1 .AND. RR .GT. STPALT) 

66 IGRAF = IGRAF + 1 

67 IF (IGRAPH .GT. 19999) THEN 

68 PRINT *, ’TOO MANY DATA POINTS’ 

69 STPALT = 1000000.0 

70 GOTO 200 

71 ENDIF 

72 GRAFX(IGRAF) = X( 1 ) 

73 GRAFY (IGRAF) = X(2) 

74 W = QSQRT(X(3)**2. + X(4)**2. ) 

75 RR =QSQRT((X(1) + DE)**2. + X(2)**2.) 

76 RMT = QSQRT((X(1) - DM)**2. + X(2)**2.) 

77 IF (X(1).LT.210000.0) THEN 

78 xf=X(l)+DE 

79 dt = INT(.01 * RR) 

80 IF (RR.LT.7800.0.AND.DIRECT.EQ.1) dt=5. 

8 1 IF( RR.GT.40000.0. AND.DIRECT.EQ.O) THEN 

82 CALL JACOBI(X,CJ,RR,EMDIST JRMT JMUN,W ,NW) 

83 ELSE 

84 CJ = 10. 

85 END IF 

86 ELSE 

87 C LUNAR PORTION 

88 xf = X(1)-DM 

89 dt = QFLOAT(INT(.02 * RMT)) 

90 IF(RMT.LT.2200.0.AND.DIRECT.EQ.0) dt=3. 

91 IF(RMT.GT.6000.0.AND.DIRECT.EQ.l) THEN 

92 CALL JACOBI(X, CJ,RR,EMDIST,RMT ) MUN,VV,NW) 

93 ELSE 

94 CJ = 10.0 
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95 END IF 

96 END IF 

97 theta = QATAN2( X(2), xf) 

98 IF (theta .LT. - (PI/2.0) ) theta = theta + 2. * PI 

99 ALPHA = theta + PI/ 2.0 

100 TT = TT + dt 

1 0 1 SCMASSV = SCMASSV - Mdot * dt * QFLOAT ( TH ) 

102 IF ( SCMASS/SCMASSV .LT. 0.0 ) THEN 

103 STPALT = 10000000.0 

104 GOTO 200 

105 ENDIF 

106 DELT_V= gravity *Isp *QLOG(SCMASS/SCMASSV) 

107 ACCEL = THRUST / SCMASSV 

108 CALL RUK4 (X, NUM, dt, MU,ACCEL1 ,DM,DE,MUM, 

109 + MUE,TH,VR,thrst,CJ,ACCEL,RANGE, DIRECT, ALPHA, 

110 * RANGEOFF, theta, RJAC,VTN,DIST,VRR,VRN,VRO, IT, VT) 

111 IF (counter.GE.5) THEN 

112 AU = 1 

113 WRTTE(AU,17)TT/ 3600.0, W,RR 

114 17 FORMAT (’OTime Hapsed = ’,F8.2,’hrs.V,’ Velocity = ’, 

115 + F9.4, ’km/s’,/,’ Dist. Earth = ’,F12.2, ’km’) 

1 16 WRITE (AU,171) RMT, SCMASS-SCMASSV, DELT_V 

117 171 FORMAT (’ Dist. Moon = ’,F12.2, ’km’,/,’ Prop, mass = ’, 

118 + F12.2,’kg’,/,’ Deltvel. = ’.F10.5, ’km/s’) 

119 WRITE (AU, 172) CJ 

120 172 FORMAT (’ JACOBIAN = ’,E15.8) 

121 counter=0 

122 ELSE 

123 counter = counter + 1 

124 ENDIF 

125 200 END DO 

126 print *, ’Orbit Completed.’ 

127 C 

128 C DRAW GRAPHICS 

129 C 

1 30 CALL GATTRI ( 1 ,0, 1 .0) 

131 CALL GATTRI (2, 0,1.0) 

132 CALL GATTRI (3, 0,1.0) 

1 33 CALL GATTRI (4,5, 1 .0) 

134 CALL GATTRI (5, 5, 1.0) 

135 CALL GATTRI (6, 5, 1.0) 

136 CALL GATTRI (7, 5,1.0) 

137 CALL GATTRI (11,5,1.0) 

138 CALL GATTRI (12,5,1.0) 

139 CALL GCHART (1,5, TITLED, 24) 

140 CALLGAXIS (1,0, -250000.0, 500000.0,0, ’Distance with respect to 

141 * bary center ’, 37, -275000.0, 325000.0,0, ’km’, 2) 
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142 RAD = 6371.23 

~ 143 DO 38 IP ASS =1,2 

144 IF (IPASS .EQ. 2) RAD = 1739.35 

145 DO 28 IN =1,200 

146 INM1 = IN - 1 

147 XPLOT(IN) = RAD * QCOS( QFLOAT (INM1) * 10. * PI/180.) 

148 YPLOT(IN) = RAD * QSIN( QFLOAT (INM1) * 10. * PI/180.) 

149 IF ( IPASS .EQ. 1 ) XPLOT ( IN ) = XPLOT ( IN ) - DE 

150 IF ( IPASS .EQ. 2 ) XPLOT ( IN ) = DM + XPLOT ( IN ) 

151 28 CONTINUE 

152 CALLJOPEN 

153 CALL JCOLOR ( 4 ) 

1 54 CALL GCURVE (XPLOT, YPLOT, 200, 0,0,0) 

-155 CALL JCLOSE 

156 38 CONTINUE 

157 CALLJOPEN 

-158 CALL JCOLOR ( 6 ) 

1 59 CALL GCURVE ( XPLOT, YPLOT, 200, 0, 0, 0 ) 

160 CALL JCOLOR ( 2 ) 

-161 CALL GCURVE ( GRAFX, GRAFY, IGRAF, 0, 0, 0 ) 

162 CALL JCLOSE 

163 CLOSE ( UNIT = 1) 

-164 CALL JPAUSE ( 1 ) 

165 CALL JDEVOF ( 1 ) 

166 CALL JDEND ( 1 ) 

-167 STOP 

168 END 
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CURVE Subroutine Description 


The CURVE subroutine is a stand alone linear, logarithmic, power, and exponential curve fitting 
routine that was taken from a Public Domain library and modified to return the class (i.e. linear, 
log., power, or exp.) of curve that best fit the input data. 

CURVE Subroutine Passed Variables 

CALL CURVE (PR(), PVT(), J, VT()) 

SUB CURVE (X(), Y(), N, VT()) 

N The number of data points. 

VT() The array of the coefficients of the curve fit. 

X() The array of the sampled radial positions. 

Y() The array of the tangential velocity component at each sampled position. 
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BASIC CODE 
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DEFDBL A-H, J-Z 

T V ?ror iNG by Don McDade - ModiIi * d by *«—y« 

REDIM A«). BH>, RU)' 

Calculate curves 

SX ■■ 0: SY — 0: SXY — 0: SXSQ ■» 0:. SYSO — 0* SXJ « 0- «;y 7 « n. cw*t n 
SXJSQ - 0 : SY JSQ - 0 :• SXX - 0 : SYK - 0 : SXYK - 0 : SXKSQ - 0 • SYKSO -■ 0 
I -'I 5 ™ • 0: SX ™ - - 0: SXMSQ - 0: 3^0: X - 0 m°- 0 

SX - SX + X (I) : SY - SY + Y(I) : SXY - SXY + X(I) * Y(I) 

X(I): SYSQ - SYSQ + Y(I) * Y(I) 'linear 
- J + 1: LY - LOG (Y (I) ) : SXJ - SXJ + X(I): SYJ - SYJ + tv . 
LY: SXJSQ - SXJSQ + X(I) 5YJ + LYl 


SXSQ - SXSQ + X(I) * 


0 

+ 


THEN J 
X(I) * 


THEN K - K + 1: LX - LOG(X(I)): 
LX * Y (I) : SXKSQ - SXKSQ + LX * 


X (I) : SYJSQ - SYJSQ + LY * LY'4xpo 

SXK - SXX + LX: SYK - SYK + Y(I) • SX 
LX: SYKSQ - SYKSQ + Y(I) * Y(I)'ioga 


AND Y (I) >0 THEN M - M + 1: SXM - SXM + LX: SYM - SYM + LY- 
LY: SXMSQ - SXMSQ + LX * LX: SYMSQ - SYMSQ + LY * LY' power 


* (N * SYSQ - SY 
B ( 1 ) ; "x"; 


IF Y (I) > 

YJ - SXYJ 
nential 
IF X(I) > 

YK - SXYK 
rithznic 
IF X(I) > 

SXYM + LX 
NEXT I 

A(l) - (SY * SXSQ - SX * SXY) / (N * SXSQ - SX * SX) 

B 1 “ (N * SXY - SX * SY) / (N * SXSQ - SX * SX) 

R(l) - (N * SXY - SX * SY) / SQR ( (N * SXSQ - SX * SX) 

PRINT SPACES (39): IF B(l) >- 0 THEN PRINT "y-"; A(l); 

T "y-"; A (1 ) ; B(l); "x"; * 1 

LOCATE CSRLIN, 56: PRINT "R-"; R(l) 

IF J < 2 THEN A (2 ) - 0: B(2) - 0: R(2) - 0 

n ?! " ^ P J ( ? YJ * SXJSQ - SXJ * SXYJ) / (J * SXJSQ - SXJ 
® i ^ - SXJ * SYJ) / (J » SXJSQ - SXJ * SXJ) 

R syj)7 J SXYJ ~ SXJ * SYJ) / SQR((J * SXJS( 2 - sxj * sxj) 

Tdy'-iVtivL if ;-; 3 ; ““ TE CSRLIN - 56: print • 

;jl| : ’ ---- : 3?s us? * “•» 

R SYM)” M SXYM ~ SXM * SYM) / SQR({M * SXMSQ - SXM * SXM) * (M 

PRINT "y-"; A (3); "x~("; B(3); ")"; : LOCATE CSRLIN, 56: PRINT -R 
THEN A ( 4 ) - 0: B(4) * 0: R ( 4 ) * 0 

(SYK * SXKSQ - SXK * SXYK) / (K * SXKSQ - SXK * SXK) 


SX 


SXYM - 


* SY)) 
ELSE PRIN 


SXJ) ) 
(J 


(K 

(K 


SXYK - SXK 
SXYK - SXK 


SYK) 

SYK) 


0 THEN PRINT "y« 


IF K < 

A (4 ) - 
B (4 ) - 
R(4) - 
SYK) ) 

IF B (4) >- 
(4); "In x ,- 

LOCATE CSRLIN, 56: PRINT 
MAXB - ABS (R ( 1 ) ) : VTN - 1 
FOR SORT — 1 TO 4 

MAX - ABS (R (SORT)) : CC - SORT 

IF MAX > MAXB THEN MAXB - MAX- VTN 
NEXT SORT 

VT(1) - A (VTN) : VT(2) - B (VTN) 

PRINT "NUMBER"; VTN; ' "WAS CHOSEN- 
END SUB 


* SYJSQ - SYJ * 
R-"; R (2) 

* SYMSQ - SYM * 

R(3) 


/ (K * SXKSQ - SXK * SXK) 

/ SQR ( (K * SXKSQ - SXK * SXK) 


A (4 ) ; 
’ R(4) 


+"; B ( 4 ) ; -In x", 


(K * SYXSQ - SYK * 
ELSE PRINT "y="; A (4); E 


CC 
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FORTRAN CODE 
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172 c SQUARES CURVE FITTING by Don McDade, Modified by David Korsmeyer 

173 SUBROUTINE CURVE (X, Y, N, VT,VR0,VTN) 

1 74 IMPLICIT REAL* 1 6 ( A-H.O-Z) 

1 75 IMPLICIT INTEGER (I-N) 

176 INTEGER VTN,AU,CC 

177 REAL* 16 MAXB,MAX,LX,LY 

178 DIMENSION A(4), B(4), R(4), VT(5), X(2100), Y(2100) 

179 AU = 1 

180 c Calculate curves 

181 SX = 0.0 

182 SY = 0.0 

183 SXY = 0.0 

184 SXSQ = 0.0 

185 SYSQ = 0.0 

186 SXJ = 0.0 

187 SYJ = 0.0 

188 SXYJ = 0.0 

189 SXJSQ = 0.0 

190 SYJSQ = 0.0 

191 SXK = 0.0 

192 SYK = 0.0 

193 SXYK = 0.0 

194 SXKSQ = 0.0 

195 SYKSQ = 0.0 

196 SXM = 0.0 

197 SYM = 0.0 

198 SXYM = 0.0 

199 SXMSQ = 0.0 

200 SYMSQ = 0.0 

201 J = 0 

202 K = 0 

203 M = 0 

204 DO 10 1 = 1, N 

205 SX = SX + X(I) 

206 SY = SY + Y(I) 

207 SXY = SXY + X(I) * Y(I) 

208 SXSQ = SXSQ + X(I) * X(I) 

209 SYSQ = SYSQ + Y(I) * Y(I) 

210 IF( Y(I).GT. 0.0) THEN 

211 J=J+ 1 

212 LY = LOG(Y(I)) 

213 SXJ = SXJ + X(I) 

214 SYJ = SYJ + LY 

215 SXYJ = SXYJ + X(I) * LY 

216 SXJSQ = SXJSQ + X(I) * X(I) 

217 SYJSQ = SYJSQ + LY * LY 

218 END IF 
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219 

220 
221 
222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 
261 
262 

263 

264 

265 


IF (X(I).GT.O.O) THEN 
K = K+ 1 
LX = LOG(X(I)) 

SXK = SXK + LX 
SYK = SYK + Y(I) 

SXYK = SXYK + LX * Y(I) 

SXKSQ = SXKSQ + LX * LX 
SYKSQ = SYKSQ + Y(I) * Y(I) 

END IF 

IF( X(I).GT.O.O. AND. Y(I)GT.O.O) THEN 
M = M + 1 
SXM = SXM + LX 
SYM = SYM + LY 
SXYM = SXYM + LX * LY 
SXMSQ = SXMSQ + LX * LX 
SYMSQ = SYMSQ + LY * LY 
END IF 

10 CONTINUE 

A(l) = (SY * SXSQ - SX * SXY) / (QFLOAT(N) * SXSQ - SX * SX) 

B(l) = (QFLOAT(N) * SXY - SX * SY) / (QFLOAT(N) * SXSQ - SX * SX) 
R(l) = (QFLOAT(N) * SXY - SX * SY) / QSQRT((QFLOAT(N) * 

+SXSQ - SX * SX) * (N * SYSQ - SY * SY)) 

IF(B(1).GE. 0.0) THEN 
WRITE(AU,17) A(l), B(l), R(l) 

ELSE 

WRITE(AU,27) A(l), B(l), R(l) 

END IF 

17 FORMAT(’OY=\E15.8,’+\E15.8,’ X\T50, ’R = \E15.8) 

27 FORMAT(’OY=’,2E15.8,’ X’,T50, ’R = ’ , E15.8) 

IF (J.LT.2) THEN 
A(2) = 0.0 
B(2) = 0.0 
R(2) = 0.0 
END IF 

A(2) = EXP((SYJ * SXJSQ - SXJ * SXYJ) / 

+(QFLOAT(J) * SXJSQ - SXJ * SXJ)) 

B(2) = (QFLOAT(J) * SXYJ - SXJ * SYJ) / 

+ (QFLOAT(J) * SXJSQ-SXJ**2.) 

R(2) = (QFLOAT(J) * SXYJ - SXJ * SYJ) / QSQRT((QFLOAT(J) * 
+SXJSQ - SXJ * SXJ) * (QFLOAT(J) * SYJSQ - SYJ * SYJ)) 
WRITE(AU,37) A(2), B(2), R(2) 

37 FORMATS y=\ E15.8, ’e (’, E15.8, ’x)\T50,’R = \E15.8) 

IF (M.LT.2) THEN 
A(3) = 0.0 
B(3) = 0.0 
R(3) = 0.0 
END IF 
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266 A(3) = EXP((SYM * SXMSQ - SXM * SXYM) /(M*SXMSQ - SXM * SXM)) 

267 B(3) = (QFLOAT(M) * SXYM - SXM * SYM) / 

268 +(QFLOAT(M) * SXMSQ - SXM * SXM) 

269 R(3) = (QFLOAT(M) * SXYM - SXM * SYM) / QS QRT((QFLOAT(M) * SXMSQ - 

270 +SXM * SXM) * (QFLOAT(M) * SYMSQ - SYM * SYM)) 

271 WRTTE(AU,47) A(3), B(3), R(3) 

272 47 FORMAT(’ y=’, E15.8, ’x (’,E15.8, ’)’,T50,’R = \E15.8 ) 

273 IF (K.LT.2) THEN 

274 A(4) = 0.0 

275 B(4) = 0.0 

276 R(4) = 0.0 

277 END IF 

278 A(4) = (SYK *SXKSQ - SXK *SXYK)/(QFLOAT(K) * SXKSQ -SXK * SXK) 

279 B(4) = (QFLOAT(K) * SXYK -SXK * SYK) / (QFLOAT(K) * SXKSQ - 

280 +SXK * SXK) 

28 1 R(4) = (QFLOAT(K) * SXYK -SXK * SYK)/QSQRT((QFLOAT(K) *SXKSQ - 

282 +SXK * SXK) * (QFLOAT(K) * SYKSQ - SYK * SYK)) 

283 IF (B(4).GE.0.0) THEN 

284 WRITE(AU,57) A(4), B(4), R(4) 

285 57 FORMAT(’ y=’,E15.8, ’+’,E15.8, ’In x’,T50, ’R = \E15.8) 

286 ELSE 

287 WRITE(AU,67) A(4), B(4), R(4) 

288 67 FORMAT(’ y=\E15.8,E15.8, ’in x’,T50,’R = ’.E15.8 ) 

289 ENDIF 

290 MAXB = QABS(R( 1 )) 

291 VTN = 1 

292 DO 20 ISORT =1,4 

293 MAX = QABS(R(ISORT)) 

294 CC = ISORT 

295 IF (MAX.GT. MAXB) THEN 

296 MAXB = MAX 

297 VTN = CC 

298 END IF 

299 20 CONTINUE 

300 VT(1) = A(VTN) 

301 VT(2) = B(VTN) 

302 WRITE(AU,1 17) VTN 

303 1 17 FORMAT (’0NUMBER ’,14/ WAS CHOSEN’) 

304 RETURN 

305 END 
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DERI Subroutine Description 


The DERI subroutine is called by the Runge-Kutta integration routine from SPIRAL. It 
contains the equations of motion for a two-body system. The only perturbation force is that of 
the spacecraft’s engine thrust. This thrust is incorporated in the equations of motion in the form 
of accelerations, a, and a,. 


x 



+ 


a 


x 


The acceleration of the spacecraft due to the thrust of the propulsion system was determined by 
calling the subroutine GUIDE. This subroutine returned the values of GDI and GD2 which are 
the magnitude of the thrusting acceleration in the x and y direction. 
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DERI Subroutine Internal and Share Variables 


SUB DERI (X(), DX(» STATIC 

DX() The array of the derivatives being integrated by the Runge-Kutta. 
gmr The gravitational force on the spacecraft at the given radial distance. 
R Radial distance of the spacecraft from the central body (km). 

r2 Square of the radius from the central body (km A 2). 

V The magnitude of the velocity of the spacecraft (km/s). 

X() State Vector of the spacecraft’s position and velocity (km and km/s). 
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DER Subroutine Flowchart 


DER 



RETURN 
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BASIC CODE 
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DEFDBL A-H, J-Z 

SOB DERI (X(), DX ( ) ) STATIC 

DX(1) - X (3) 

DX(2) - X (4) 

r2 - X(l) ^ 2 + X (2) - 2 
R - SQR(r2) 
gmr - -MU / r2 

V - SQR (X (3) * 2 + X (4 ) - 2) 
DX(3) - gmr * X(l) / R + ACCEL1 
DX (4 ) - gmr * X(2) / R + ACCEL1 
END SUB 


X (3) / V 
X(4) / V 
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FORTRAN CODE 
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335 SUBROUTINE DER 1 (X, DX,MU, ACCELI ) 

336 IMPLICIT REAL* 1 6 (A-Z) 

337 DIMENSION X(4),DX(4) 

338 DX(I) = X(3) 

339 DX(2) = X(4) 

340 r2 = X( 1 ) ** 2 + X(2) ** 2 

341 R = QSQRT(r2) 

342 gmr = -MU /r2 

343 V = QSQRT(X(3) ** 2 + X(4) ** 2) 

344 DX(3) = gmr * X(l) /R +ACCEL1 *X(3)/ V 

345 DX(4) = gmr * X(2) / R + ACCELI * X(4) / V 

346 RETURN 

347 END 


61 



DER Subroutine Description 


The DER subroutine contains the equations of motion for the 
body system. These are: 

x - ny - nx = -- I + 

r. r m 


y + 2nx - n 2 y 



spacecraft in the restricted three- 


] 

] 


Where a, is the guidance acceleration in the y-direction, a, is the guidance acceleration in the x- 
direction, |i is the gravitational parameter of the target planet, and r is the radial distance of the 
spacecraft from the target planet. 
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DER Subroutine Internal and Shared Variables 

SUB DER (X(), DX()) STATIC 

DX() The array of the derivatives being integrated by the Runge-Kutta. 
dx3 The second derivative of the x-component (km/s A 2). 

dx4 The second derivative of the y-component (km/s A 2). 

GDI The x-acceleration of the spacecraft from GUIDE (km/s A 2). 

GD2 The y-acceleration of the spacecraft from GUIDE (km/s A 2). 

muroe The gravitational parameter of the Earth divided by ROE. 

murom The gravitational parameter of the Moon divided by ROM. 

ROE The cubed distance of the spacecraft from the Earth in the rotating x,y coordinates 
(km A 3). 

X() State Vector of the spacecraft’s position and velocity in the rotating x,y coordinates 
(km and km/s). 

xde The x-coordinate of the Earth in the rotating x,y coordinates (km), 

xdm The x-coordinate of the Moon in the rotating x,y coordinates (km). 
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BASIC CODE 
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DEFDBL A-H, J-Z 


SUB DER (X<) 
DX(1) - X (3) 


DX ( ) ) STATIC 


xdm - X (1) - DM 


DX(2) 
xde - 
X22 - 
ROE - 
ROM - 
muroe 
dx3 - 
) 

dx4 - 

) 

CALL GUIDE (GDI, GD2, 0C()) 

DX(3) - dx3 + GDI: DX(4) - dx4 + 
END SUB 


'xdot (tan/ 3 ) 
'ydot (tan/ a) 


- X (4 ) 

X(l) + DE 
X (2) ^ 2 

(xde * xde + X22) * 1.5# 

(xdm * xdm + X22) "'1.5# 

- MUE / ROE : murom - MUM / ROM * mu/radius ratios 
-muroe * xde - murom * xdm + 2# * X(4) * NW + NW2 


' coordinates 
' y-squared 
' distance of 
' distance of 


of the Earth and Moon 

s/c from Earth (tan) 
s/c from Moon (tan) 


X(l) ' x-dbldct 


-X(2) * (muroe + murom) - 2# * X(3) * NW + NW2 * X (2) 'y-dbldot 


GD2 


(tan/s"I 

(tan/s"2 
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J08 SUBROUTINE DER (X, DX, DM, DE ,MUM , MUE ,TH ,VR, thrst, 

309 * CJ, ACCEL, RANGE, DIRECT, ALPHA, RANGEOFF.THETA, 

310 * RJAC,VTN,VT,DIST,VRR,VRN,VRO,TT) 

1 1 1 IMPLICIT REAL* 16 (A-H.J-Z) 

312 DIMENSION DX(4),X(4),VR(10),VT(5), RJAC(3) 

3 1 3 INTEGER VTN,VRN,DIRECT,RANGEOFF,TH, thrst 

114 NW= 2.6653 14572E-06 

315 nw2 = nw**2. 

316 DX(1) = X(3) 

117 DX(2) = X(4) 

318 xde = X(l) + DE 

319 xdm = X(l) - DM 

120 X22 = X(2) ** 2. 

321 ROE = (xde * xde + X22) ** 1.5 

322 ROM = (xdm * xdm + X22) ** 1 .5 

•23 muroe = MUE / ROE 

324 murom = MUM / ROM 

325 dx3 = -muroe * xde - murom *xdm +2 * X(4) * NW + NW2 * X(l) 

26 dx4 = -X(2) * (muroe + murom) - 2 * X(3) * NW + NW2 * X(2) 

327 CALL GUIDE(GD 1 , GD2, X, DIRECT, PI,ACCEL, DM, ALPHA.DE.RANGEOFF, 

328 + RANGE, thrst, CJ, theta, RJAC,TH,VTN,VT,DIST,VR,VRR, 

29 + VRN,VR0,TT) 

-^30 DX(3) = dx3 + GDI 

331 DX(4) = dx4 + GD2 

32 RETURN 

-o33 END 
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ERASE Subroutine Description 


This subroutine is used to reinitialize the velocity component arrays. The arrays are used to 
provide curve fit data to the POLYFIT subroutine. This subroutine does not exist in the BASIC 
version of CISGRAPH. 


ERASE Subroutine Variables 


VRO First point of the radial velocity array used to produce guidance curves. 
VR() Radial velocity array used to produce guidance curves. 

VT() Tangential velocity array used to produce guidance curves. 
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FORTRAN CODE 
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349 SUBROUTINE ERASE(VR, VT, VRO) 

350 REAL* 16 VR0,VR(10),VT(5) 

351 VR0= 0.0 

352 DO 10 1=1,5 

353 VR(I) = 0.0 

354 VR(I+5) = 0.0 

355 VT(I) = 0.0 

356 10 CONTINUE 

357 RETURN 

358 END 



GUIDE Subroutine Description 


The GUIDE subroutine controls the direction and magnitude of the spacecraft’s thrusting 
acceleration. A set of guidance parameters are determined depending on the direction of the 
spacecraft’s trajectory. When the spacecraft is in the departure phase of the trajectory the 
acceleration is along the velocity vector of the spacecraft. When the spacecraft Jacobian value 
passes the first control Jacobian, Jacl, the program checks the quadrant the spacecraft is in, and 
leaves the thrust on, if it is in the appropriate quadrant; turns it off if not. If the spacecraft’s 
Jacobian value has passed the second control Jacobian, Jac2, then the spacecraft returns to 
continuously thrusting along the velocity vector. When the spacecraft’s Jacobian has passed the 
third control Jacobian, Jac3, the spacecraft’s thrust is turned off and the spacecraft drifts until its 
distance from the departure planet has passed the value of RANGE. This initiates the capture 
phase of the trajectory. The tangential and radial components of the reference trajectory are 
calculated from their parametric functions and the direction and magnitude of the capture 
acceleration is determined. 
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GUIDE Subroutine Internal and Shared Variables 

SUB GUIDE (GDI, GD2, X()) 

ACCEL Acceleration of the spacecraft during the cislunar trajectory (km/s A 2). 

ACXT The tangential acceleration of the spacecraft in the x -direction (km/s A 2). 

ACYT The tangential acceleration of the spacecraft in the y-direction (km/s A 2). 

alph The angle between the velocity vector and the tangential component of velocity 

(radians). 

ALPHA theta + PI/2 

am Magnitude of the capture guidance acceleration (km/s A 2). 

ami The capture guidance acceleration in the x-direction (km/s A 2). 

amy The capture guidance acceleration in the y-direction (km/s A 2). 

ANGLE 1 The angle defining the beginning of the third quadrant area around the departure 
planet during the departure phase of the trajectory generation (radians). 

ANGLE2 The angle defining the end of the third quadrant area around the departure planet 
during the departure phase of the trajectory generation (radians). 

CJ Instantaneous Jacobian constant of the spacecraft. 

DIRECTION Numeric indicator of the direction of the trajectory generation, Earth to Moon (0) 
or Moon to Earth (1). 

DIST The distance the spacecraft is away from the capture planet (km). 

dt Integration step size (seconds). 

GDI The x- acceleration of the spacecraft from the guidance control equations 

(km/s A 2). 

GD2 The y-acceleration of the spacecraft from the guidance control equations 

(km/s A 2). 

rl The distance from the capture planet (km/1000). 

RANGE Range from initial planet that the capture phase is initiated (km). 
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RANGEOFF Flag for GUIDE subprogram indicating if the spacecraft has entered the capture 
phase of the trajectory generation. 

RJAC() The array of the Jacobian constants used as controls for the departure portion of 
the trajectory generation. 


tc 

TH 

theta 

thrst 

TT 

VMAG 

VR() 

VRN 

VRR 

VT() 

VTN 

VTT 

VTV 

X() 

XACCEL 

XRANGE 


Empirical time constant used to change the control velocity difference into an 
acceleration (seconds). 

Numeric indicator of the spacecraft’s thrust, on (1) or off (0). 

Angle between the radius vector to the spacecraft from the controlling gravita- 
tional body and the x-axis. Dependent upon xf. 

Flag for the spacecraft on its spiral escape indicating the passage of the second 
control Jacobian, (0) off (1) on. 

Total time of trajectory generation (seconds). 

Magnitude of the velocity of the spacecraft (km/s A 2). 

The array of parameterized radial velocity component coefficients. 

Degree of the radial velocity polynomial curve fit. 

The parameterized radial velocity magnitude (km/s). 

The array of the parameterized tangential velocity component coefficients. 

Class of equation for tangential velocity parameterization, linear (1), exponential 
(2), power (3), or logarithmic (4). 

The parameterized tangential velocity magnitude (km/s). 

Loop variable for the radial velocity polynomial. 

State Vector of the spacecraft’s position and velocity in the rotating x,y coor- 
dinates (km and km/s). 

The magnitude (+ or -) of the ACCEL used (km/s A 2). 

The magnitude of the distance in the x - direction the spacecraft is from the 
departure planet (km). 
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GUIDE Subroutine Flowchart 
GUIDE 



RETURN 
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BASIC CODE 
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DEFDBL A-H, J-Z 

SUB GUIDE (GDI, GD2, X()) 'Guidance subprogram 

SHARED TT, theta, ALPHA, CJ, dt, SQ4ASSV, thrst 
SHARED VT () , VR(), DIRECTION, ACCEL, RJAC(), RANGEOFF 
SELECT CASE DIRECTION 

CASE 0 » earth to moon 


ANGLE 1 - PI: ANGLE 2 - 1.5# 
DIST - SQR ( (X (1) - DM) ~ 2 
CASE 1 

ANGLE 1 - -6#: ANGLE2 - 0#: 
DIST - SQR ( (X (1) + DE) ~ 2 
ALPHA + PI 
END SELECT 

VMAG - SQR (X (3) A 2 + X(4) A 2) 
ACXT - ACCEL * X (3 ) / VMAG: ACYT - 
gential 

IF RANGEOFF - 1 THEN RANGE - 50000 
IF XRANGE < RANGE THEN 


* PI: XACCEL - ACCEL 

+ X(2) - 2): XRANGE - X(l): alph - ALPHA 
'moon to earth 
XACCEL - -ACCEL 

+ X(2) ^ 2): XRANGE - ABS (DM - X(l)) : alph 


ACCEL * X(4) / VMAG 'components of accel, T 


IF thrst - 0 THEN 

IF CJ < RJAC(l) THEN 

IF theta > ANGLE1 AND theta < ANGLE2 THEN 


IF CJ < RJAC (2 ) THEN thrst - 1 
TH - 1# 


ELSE 


TH - 0# 

END IF 

ELSE 
TH - 1# 

END IF 

ELSEIF CJ > RJAC (3) THEN 
TH - 1# 

ELSE 
TH - 0# 

END IF 

GDI - ACXT * TH: GD2 — ACYT * TH 'thrusting tangentially spiral out 


IF RANGEOFF - 0 THEN LINE (X(l), X(2) + 5000)- (X(l), X(2) - 5000), 
RANGEOFF - 1 ' 

tc - 100#: rl- - DIST / 1000# 

SELECT CASE VTN 
CASE 1 

VTT - VT(1) + VT (2) * rl 

CASE 2 

VTT - VT (1) * EXP (VT (2) ) 

CASE 3 

VTT - VT (1) * (rl) * VT (2) 

CASE 4 

VTT - VT (1 ) + VT (2 ) * LOG (rl) 

END SELECT 
VRR - VR ( 0 ) 

FOR VTV - 1 TO VRN 

VRR - VRR + VR (VTV) * rl * VTV 

NEXT VTV 

amx - -ACXT + ((VRR * COS (theta) - VTT * COS (alph) ) - X (3 ) ) / tc 

amy - -ACYT + ((VTT * SIN(alph) + VRR * SIN (theta)) - X(4)) / tc 

am - SQR ( amx * 2 + amy ~ 2 ) 

GDI — ACCEL / am * amx 

GD2 - ACCEL / am * amy: TH - 1# 

END IF 
END SOB 


15 
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FORTRAN CODE 
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361 

362 

363 

364 

365 

366 

367 

368 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 

380 

381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 

394 

395 

396 

397 

398 

399 

400 

401 

402 

403 

404 

405 

406 

407 


SUBROUTINE GUIDE (GDI, GD2, X, DIRECT, PI, ACCEL, DM,ALPHA,DE, 
+ RANGEOFF, RANGE, thrst, CJ, theta, RJAC.TH, 

+ VTN,VT,DIST,VR,VRR,VRN,VRO,TT) 

IMPLICIT REAL* 16 (A-H,J-Z) 

IMPLICIT INTEGER (I) 

INTEGER VTN,VRN, DIRECT, RANGEOFF, TH, thrst 
DIMENSION X(4),RJAC(3),VT(5),VR(10) 

C IF (TT/3600. .GT. 312. )PRINT *, ’VRO ’, VR0 
PI = 3.1415926535 
IF (DIRECT.EQ.O) THEN 
ANGLE1 = PI 
ANGLE2 = 1.5 * PI 
XACCEL = ACCEL 

DIST = QSQRT((X(1) - DM) ** 2. + X(2) ** 2.) 

XRANGE = X(l) 
alph = ALPHA 
ELSE 

ANGLE1 = -6.0 
ANGLE2 = 0.0 
XACCEL = -ACCEL 

DIST = QSQRT((X(1) + DE) ** 2. + X(2) ** 2.) 

XRANGE = QABS(DM-X(1)) 

Alph = ALPHA + PI 
ENDIF 

VMAG = QSQRT(X(3) ** 2. + X(4) ** 2.) 

ACXT = ACCEL * X(3) / VMAG 
ACYT = ACCEL * X(4) / VMAG 
IF (RANGEOFF.EQ.l) RANGE = 50000 
IF (XRANGE.LT.RANGE) THEN 
IF ( thrst .EQ. 0 ) THEN 
IF ( a .LT. RJAC(l) ) THEN 

IF ( theta .GT. ANGLE 1 .AND. theta .LT. ANGLE2 ) THEN 
IF ( CJ. LT. RJAC(2) ) thrst = 1 
TH=1 
ELSE 
TH=0 
ENDIF 
ELSE 
TH=1 
ENDIF 

ELSE IF (CJ.GT.RJAC(3)) THEN 
TH=1 
ELSE 
TH=0 
ENDIF 

GD1=ACXT * QFLOAT(TH) 

GD2 = ACYT * QFLOAT(TH) 
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408 ELSE 

-409 RANGEOFF=l 

410 tc = 100. 

411 rl = DIST / 1000. 

412 IF (VTN.EQ. 1 ) VTT=VT( 1 ) +VT(2) * rl 

413 IF (VTN.EQ. 2) VTT=VT(1) * EXP(VT(2)) 

414 IF (VTN.EQ.3) VTT = VT(1) * rl **VT(2) 

415 VRR = VRO 

416 IF (VTN.EQ.4) VTT = VT( 1 ) + VT(2) * LOG(rl) 

417 DO 10 VTV= 1,VRN 

418 VRR= VRR + VR(VTV) * rl ** VTV 

419 10 CONTINUE 

420 amx = -ACXT + ((VRR * QCOS(theta) - VTT*QCOS(alph))-X(3))/tc 

421 amy = -acyt + ((VTT*QSIN(alph) + VRR*QSIN(theta))-X(4))/tc 

422 am=QSQRT(amx**2. + amy**2.) 

423 GDI = ACCEL/ am * amx 

_ 424 GD2 = ACCEL/ am * amy 

425 TH = 1 

426 ENDIF 

_ 427 RETURN 

428 END 
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IO Subroutine Description 


The IO subroutine handles the initial input of spacecraft characteristics, control variables, and 
various use choices. The initial output screen formatting and graphics setup is also performed in 
this subroutine. The spacecraft’s initial operating characteristics, such as initial mass, and the 
mass flow rate of the propulsion system, are selected in addition to the direction of the trajectory 
generation. Then the user is asked whether a new set of guidance velocity parametrics should be 
generated based upon the input spacecraft characteristics. If the response is yes, the program 
control is passed to another subroutine called SPIRAL. This subroutine generates the parametric 
velocity profiles. Upon completion of SPIRAL or if the response to generating velocity 
parametrics is no, the IO subroutine prompts the user to input the spacecraft’s initial position and 
velocity. The next question asks if the user would like to modify the control Jacobians and 
Range. These control values are used by the program to control when the spacecraft escapes 
from its outbound spiral and begins the capture phase of the trajectory. The default values are 
presented and the user can modify any of them. After all of the inputs and responses have to be 
recorded the IO subroutines sets up the screen for the trajectory generation and returns control to 
the Main program. 
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IO Subroutine Internal and Shared Variables 


SUB IO (VR(), X(), DIRECTION, RJAC()) STATIC 

DIRECTION Numeric indicator of the direction of the trajectory generation, Earth to Moon (0) 
or Moon to Earth (1). 

rl Initial radius value, input from keyboard (km). 

RJAC() The array of the Jacobian constants used as controls for the departure portion of 
the trajectory generation. 

thet Angle, in radians, of theta. 

theta Initial angle, in degrees, spacecraft is from - x-axis, input from keyboard. 

title$ String Variable for the tide of the output screen. 

vexp Velocity magnitude of spacecraft, input from keyboard (km). 

VR() The array of the parameterized radial velocity component coefficients. 

VT() The array of the parameterized tangential velocity component coefficients. 

X() State Vector of the spacecraft's position and velocity in the rotating x,y coordi- 

nates (km and km/s). 
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IO Subroutine Flowchart 
10 



Return 
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BASIC CODE 
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DEFDBL A-H, K-Z 

SUB 10 (VR<), VT(), X(), DIRECTION, RJACO) STATIC 
CLS 

SCREEN 0 

PRINT m ZDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDL 
DDD?~ 

PRINT "J Trajectory Generation Model for Low-Thrust OTVs in Cislunar Space 
3* 

PRINT "J Using a Thrusting Control Algorithm in the Restricted three-body 

3 ” 

PRINT " 3 formulation of the Earth-Moon system. {Gus Babb copy} - DJK, 4/l/8 c 
J" 


PRINT ” SDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD- 
DDDY* 

DBOX 6, 1, 11, 80 

LOCATE 7, 3: PRINT "Spacecraft intial mass SCMASS; 

LOCATE 8, 3: PRINT "Specific Impulse of engine -", Isp; 

LOCATE 9, 3: PRINT "Mass flow rate of engines - ", Mdot 

LOCATE 10, 3: INPUT "Do you wish to specify s/c characteristics? (y or n)", B* 
IF UCASE$(B$) - "Y" THEN 
DBOX 11, 1, 16, 80 

LOCATE 12, 3: INPUT "Spacecraft intial mass -", SCMASS 

LOCATE 13, 3: INPUT "Specific Impulse of engine -", Isp 

LOCATE 14, 3: INPUT "Mass flow rate of engines - ", Mdot 

LOCATE 15, 3: INPUT "Degree of polynomial curve fit (2-7) ", VRN 

END IF 

THRUST - gravity * Isp * Mdot 'thrust dependent on mass flow and isp 

PRINT 

INPUT "Starting Orbit about the Earth or the Moon? (e or m) ", A$ 

IF UCASES (AS) - "M" THEN 

DIRECTION - 1 'flag indicates s/c going moon to earth 

RANGE - 110000 

titleS - "Moon to Earth Trajectory" 

INPUT "Do you want to generate parmetric velocity curves? ”, AS 
IF UCASES(AS) - "Y" THEN 
ERASE VT, VR, X 

X(0) - 0#: X (2 ) - 7178#: X<3) SQR (MUE / X(2)): X<4) - 0# 

SPIRAL X(), DIRECTION, 225000, VT(), VR() 


ELSE 

VR (0) - 10.92448764296168# 

VR (1 ) - -.6486459794313275# 

VR (2)- - .0233450835020553# 

VR(3) 4 . 502555723665775D-04 

VR (4) - 4 . 767029852095007D-06 
VR (5 ) - -2 . 760604361514407D-08 
VR (6) - 8 . 185026033125225D-11 
VR (7 ) - -9 . 6764120959764D-14 
VT (1) - 3. 880360674664431D-02 

VT (2) 4 . 680399272731233D-03 

VTN - 1 


END IF 

INPUT "Input the radius from the Moon's center ", rl 
INPUT "Input angle (deg) from -x axis ", theta 
INPUT "Do you wish to specify the velocity? (y or n) ", B$ 
IF UCASES(BS) - "Y" THEN 

INPUT "Input absolute velocity ", vexp 

ELSE 


vexp - -SQR (MUM / rl) 

END IF 

thet - theta * PI / 180# 

X(l) - rl * COS (thet) + DM 
X (2) - rl * SIN (thet) 

X (3) - -vexp * SIN (thet) 

X (4) - vexp * COS (thet) 

INPUT "Modify Control Jacobians and Range? ", C$ 
IF UCASES (CS) - "Y" THEN 
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PRINT "Default Values 3.05, 2.82, 2.50" 

PRINT " Jacl “ ", RJAC(l): INPUT RJAC(l) 

PRINT " Jac2 - ", RJAC ( 2 ) : INPUT RJAC(2) 

PRINT "Jac3 - ", R JAC ( 3 ) : INPUT RJAC(3) 

PRINT "Range - ", RANGE: INPUT RANGE 

ELSE 

RJAC(l) - 3.05#: RJAC (2) - 2.82#: RJAC(3) - 2.5# 

END IF 

ELSE 

DIRECTION - 0 'flag indicates s/c going earth to moon 

RANGE - 255000 

titles “ "Earth to Moon Trajectory" 

INPUT "Do you want to generate parmetric velocity curves? ", A$ 

IF UCASE$(A$) - "Y" THEN 
ERASE VT, VR, X 

X(0) - 0#: X (2) - 1838#: X(3) - -SQR (MUM / X (2 ) ) : X ( 4 } - 0# 
SPIRAL X(), DIRECTION, 120000, VT(), VR() 

ELSE 

VR (0 ) - 2.083801929462969* 

VR (1 ) - -.328978346551485# 

VR (2) - 2. 759195072751214D-02 
VR (3) - -1 . 13520698564 688D-03 
VR ( 4 ) - 2 . 4 54 2 9 61 3 154 4 3 68D- 05 
VR (5) - -2.850460877535156D-07 
VR (6) - 1. 681222735220663D-09 
VR (7) - -3. 948910484014757D-12 
VT (1 ) - 2.2187* 

VT (2 ) - -.5012# 

VTN - 3 

END IF 

PRINT "Input the radius from the Earth's center " 

INPUT rl : PRINT "Input angle (deg) from -x axis " 

INPUT theta: INPUT "Do you wish to specify the velocity? (y or n) ", 3$ 
IF UCASES(BS) - "Y" THEN 

INPUT "Input absolute velocity ", vexp 

ELSE 

vexp - SQR (MUE / rl) 

END IF 

thet - theta * PI / 180# 

X(l) = rl * COS (thet) - DE 
X (2 ) - rl * SIN (thet) 

X(3) - -vexp * SIN (thet) 

X ( 4 ) «= vexp * COS (thet) 

INPUT "Modify Control Jacobians and Range? ", C$ 

IF UCASES(CS) - "Y" THEN 

PRINT "Default Values 4.93, 4.10, 2.70" 

PRINT "Jacl - ", RJAC ( 1 ) : INPUT RJAC(l) 

PRINT " Jac2 - ", RJAC ( 2 ) : INPUT RJAC(2) 

PRINT " Jac3 - ", RJAC(3): INPUT RJAC(3) 

PRINT "Range - ", RANGE: INPUT RANGE 

ELSE 

RJAC ( 1 ) - 4.93#: RJAC(2) - 4.1#: RJAC(3) - 2.7# 

END IF 

END IF 
CLS 1 

SCREEN 2: WINDOW (-250000, -275000) - (500000, 325000) 

CIRCLE (-DE, 0), 6378 
PAINT (-DE, 0), 9, 15 
CIRCLE (DM, 0), 1734 
PAINT (DM, 0), 8, 15 
DBOX 1, 1, 4, 80 
DBOX 1, 1, 25, 80 
LOCATE 1, 28: PRINT titleS 


LOCATE 2, 2: PRINT "I nit. Mass 
kg" 


: PRINT USING "######.##"; SCMASS; : PRINT " 


LOCATE 2, 27: PRINT "Prop. Mass - "; 
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LOCATE 3, 53: PRINT "Dist. Moon - : LOCATE 2, 53: PRINT "Dist. Earth - -■ 

LOCATE 3, 2: PRINT -Time elapsed - •; 

LOCATE 3, 30: PRINT ?Vel - 
LOCATE 5, 2: PRINT -Jacobian - 

LOCATE 25, 21: PRINT "P - Pause G - Go R - Restart Q-quit"; 

END SUB 
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FORTRAN CODE 
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43 1 SUBROUTINE 10 (SCMASS,Isp,Mdot,VRN,THRUST,X, 

432 * gravity .TITLED, RANGE, DIRECT ,MUE,VR,rl, theta, RJAC, PI, 

433 * VT,VRR>lUACCELl,MUMNUM,DM,DE,dt,STPALT,TT,VTN,VRO) 

434 IMPLICIT REAL* 1 6 ( A-Z) 

435 CHARACTER* 24 TITLED 

436 character*4 AD 

437 CHARACTER* 1 BD, CD 

438 INTEGER VTN, DIRECT, NUM.VRN, AUJ 

439 DIMENSION X(4),VR(10),VT(5),RJAC(3) 

440 AU = 5 

441 WRITE (AU, 71) 

442 WRITE (AU,72) 

443 AU = 1 

444 WRITE (AU,71) 

445 WRITE (AU.72) 

446 7 1 FORMAT ( ’ ITrajectory Generation Model for Low -Thrust OTV’ ’s V* 

447 * ’ in Cislunar Space using a Thrusting Control’ J, 

448 *’ Algorithm in the Restricted three-body’) 

449 72 FORMAT (’ formulation of the Earth-Moon system. V» 

450 *’ Eagle Engineering, Inc. (LSPI - djk)’) 

451 AU = 5 

452 PRINT *, ’ ’ 

453 PRINT *, ’OUTPUT WILL GO TO FILE "CISLUNAR.OUT" ’ 

454 WRITE(AU,7) SCMASS 

455 7 FORMAT(’OSpacecraft initial mass =’,F10.2) 

456 WRTTE(AU,17) Isp 

457 17 FORMAT( ’ Specific Impulse of engine = ’, F8.2) 

458 WRTTE(AU,27) Mdot 

459 27 FORMAT (’ Mass flow rate of engines = ’,E15.8) 

460 WRTTE(AU,37) 

461 37 FORMAT (’ODo you wish to specify s/c characteristics? (y or n)’) 

462 READ *, BD 

463 IF (BD.EQ.’y’.OR.BD.EQ.’Y’) THEN 

464 PRINT *, ’ Input Spacecraft initial mass.’ 

465 READ *, SCMASS 

466 PRINT *, ’ Input Specific Impulse of engine.’ 

467 READ *,Isp 

468 PRINT *, ’ Input Mass flow rate of engines.’ 

469 READ *, Mdot 

470 PRINT *, ’ Input Degree of polynomial curve fit (2-7)’ 

471 READ *, VRN 

472 ENDIF 

473 THRUST = gravity * Isp * Mdot 

474 WRITE(AU,47) 

475 47 FORMAT (’OStarting Orbit About the Earth or Moon?’ ) 

476 READ *, AD 

477 PRINT *, ’Input Destination Altitude at which to stop processing.’ 
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478 

479 

480 

481 

482 

483 

484 

485 

486 

487 

488 

489 

490 

491 

492 

493 

494 

495 

496 

497 

498 

499 

500 

501 

502 

503 

504 

505 

506 

507 

508 

509 

510 

511 

512 

513 

514 

515 

516 

517 

518 

519 

520 

521 

522 

523 

524 


READ *, S TP ALT 
PRINT*,’ ’ 

PRINT *,’ ’ 

PRINT *,’ PLEASE WAIT, GENERATING PARAMETRIC VELOCITY EQN.’ 
IF (AD.EQ.’m’.OR.AD.EQ.’M’.OR.AD.EQ.’moon’) AD=’MOON’ 

IF (AD.EQ.’MOON’) THEN 
DIRECT = 1 
RANGE = 110000. 

TrTLED=’Moon to Earth Trajectory’ 

CALL ERASE (VR, VT, VR0) 

X(1) = 0.0 
X(2) = 7178. 

X(3)=-QSQRT(MUE/X(2» 

X(4)= 0.0 

RANGESC = 225000.0 

CALL SPIRAL (X,DIRECT,RANGESC ,VT,VR,MUE,MUM,SCMASS, 

+ dt, THRUST, TT,MdotMU,ACCELl ,NUM,DM,DE,VRN,VTN, VR0) 

PRINT *, ’Input the radius from the Moon”s Center’ 

READ *, rl 

PRINT *, ’Input angle (deg) from -x axis’ 

READ *, theta 

PRINT *, ’Do you wish to specify the velocity? (y or n)’ 

READ *, BD 

IF (BD.EQ.’y’.OR.BD.EQ.’Y’) THEN 
WRITE(AU,57) 

57 FORMAT (’OInput absolute velocity ’ ) 

READ *,VEXP 
ELSE 

vexp = -QSQRT(MUM/rl) 

ENDIF 

thet = theta * PI / 180. 

X(1 ) = rl * QCOS(thet) + DM 
X(2) = rl * QSIN(thet) 

X(3) = -vexp * QSIN(thet) 

X(4) = vexp * QCOS(thet) 

PRINT *, ’Modify Control Jacobians and Range?’ 

READ *, CD 

IF (CD.EQ.’Y’.OR.CD.EQ.’y’) THEN 
WRITE( AU,67 ) 

67 FORMAT (’ODefault Values 3.05, 2.82, 2.50’) 

PRINT *,’ JAC(l) = ’ 

READ *, RJAC(l) 

PRINT *,’ Jac2 = ’ 

READ *, RJAC(2) 

PRINT *,’ Jac3 = ’ 

READ *, RJAC(3) 

PRINT *,’ Range = ’ 
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525 

526 

527 

528 

529 

530 

531 

532 

533 

534 

535 

536 

537 

538 

539 

540 

541 

542 

543 

544 

545 

546 

547 

548 

549 

550 

551 

552 

553 

554 

555 

556 

557 

558 

559 

560 

561 

562 

563 

564 

565 

566 

567 

568 

569 

570 

571 


READ *, RANGE 
ELSE 

RJAC(1)=3.05 
RJAC(2) = 2.82 
RJAC(3) = 2.5 
ENDIF 
ELSE 

DIRECT = 0 

RANGE = 255000.0 

TITLED = ’Earth to Moon Trajectory’ 

CALL ERASE (VR,VT,VR0) 

X(1) = 0.0 
X(2) = 1838. 

X(3) = -QSQRT(MUM / X(2)) 

X(4) = 0. 

RANGESC = 120000.0 

CALL SPIRAL (X, DIRECT, RANGESC, VT, VR, MUE, MUM, 
+ SCMASS, dt, THRUST, TT, Mdot, 

+ MU,ACCEL1 ,NUM, DM , DE,VRN ,VTN,VR0) 

PRINT *, ’Input the radius from the Earth’ ’s center ’ 

READ *, rl 

PRINT *, ’Input angle (deg) from -x axis ’ 

READ *, theta 

PRINT *, ’Do you wish to specify the velocity? (y or n) ’ 

READ *, BD 

IF (BD.EQ.’Y’.OR.BD.EQ.’y’) THEN 
WRTTE(AU,77) 

77 FORMAT (’OInput absolute velocity ’ ) 

READ *, vexp 
ELSE 

vexp = QSQRT(MUE / rl) 

ENDIF 

thet = theta * PI / 180. 

X(l) = rl * QCOS(thet) -DE 
X(2) = rl * QSIN(thet) 

X(3) = -vexp * QSIN(thet) 

X(4) = vexp * QCOS(thet) 

PRINT *, ’Modify Control Jacobians and Range? ’ 

READ *, CD 

IF (CD.EQ.’Y’.OR.CD.EQ.’y’) THEN 
WRITE(AU,87) 

87 FORMAT (’ODefault Values 4.93, 4.10, 2.70’ ) 

PRINT *,’ Jacl = ’ 

READ *, RJAC(l) 

PRINT *, ’ Jac2 = ’ 

READ *, RJAC(2) 

PRINT *, ’ Jac3 = ’ 
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573 

574 

575 

576 

577 

578 

579 

580 

581 

582 


READ *,RJAC(3) 
PRINT *, ’ Range = ’ 
READ *, RANGE 
ELSE 

RJAC(l) = 4.93 
RJAC(2) = 4.1 
RJAC(3) = 2.7 
ENDIF 
ENDIF 
RETURN 
END 



JACOBI Subroutine Description 


The JACOBI subroutine calculates the instantaneous Jacobian constant of the spacecraft during 
its flight. This constant is then compared to a series of user defined control values to determine 
the necessary guidance control action. 
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JACOBI Internal and Shared Variables 


CALL JACOBI (XQ, CJ) STATIC 

CJ The instantaneous Jacobian constant of the spacecraft. 

EN The non-dimensionalized energy of the spacecraft in the three-body system. 

RMT Radial distance of the spacecraft from the center of the Moon (km). 

ROEN The non-dimensionalized distance the spacecraft is away from the Earth. 

ROMN The non-dimensionalized distance the spacecraft is away from the Moon. 

RR Radial distance of the spacecraft from the center of the Earth (km). 

VELN The non-dimensionalized velocity of the spacecraft. 

W Magnitude of the spacecraft’s velocity (km/s). 

X() State Vector of the spacecraft’s position and velocity in the rotating x,y coordinates 
(km and km/s). 

XN The non-dimensionalized x-component of the spacecraft’s position. 

YN The non-dimensionalized y-component of the spacecraft’s position. 
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BASIC CODE 


94 



DEFDBL A-H, J-Z 

SUB JACOBI (X () , CJ) STATIC 

SHARED W, RMT, RR 

ROEN - RR / EMDIST 

ROMN - RMT / EMDIST 

XN - X(l) / EMDIST 

YN - X (2) / EMDIST 

EN-XN~2 + YN-2 + 2#* (1#- MUN) / ROEN + 2# * MUN / ROMN 
VELN - W / (NW * EMDIST) 

CJ - EN - VELN * VELN + MUN * (1# - MUN) 

END SUB 
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FORTRAN CODE 
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584 SUBROUTINE JACOBI (X,CJ, RR, EMDIST, RMT, MUN, W, NW) 

585 IMPLICIT REAL* 1 6 (A-Z) 

586 DIMENSION X(4), VR(10), RJAC(3) , VT(5) 

587 ROEN = RR / EMDIST 

588 ROMN = RMT / EMDIST 

589 XN = X(l) / EMDIST 

590 YN = X(2) / EMDIST 

591 EN = XN ** 2 + YN ** 2 + 2. * (1 . - MUN) /ROEN +2.* MUN / ROMN 

592 YELN = W / (NW * EMDIST) 

593 CJ = EN - VELN * VELN + MUN * (1. - MUN) 

594 RETURN 

595 END 
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POLYFIT SUBROUTINE DESCRIPTION 


The POLYFIT subroutine is a stand along polynomial curve fitting routine that was taken from a 
Public Domain library of programs. It can fit up to seventh-order polynomial curves to the input 
data. 

POLYFIT SUBROUTINE PASSED VARIABLES 

CALL POLYFIT (PR(), PRV(), VRN, J, VR()) 

SUB POLYFIT (PR(), PV(), DEGREE, N, VR()) 

DEGREE The degree of the polynomial curve for the data. 

N The number of data points. 

PR() The array of the sampled radial positions. 

PV() The array of the radial velocity component at each sampled position. 

VR() The output array of the coefficients for the polynomial curve. 
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BASIC CODE 
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DEFDBL A-H, J-Z 

SUB POLYFIT (PRO. PV(), DEGREE, N, VR()) 

REDIM A(21), R (13, 14), t(14) 

D - DEGREE: A(l) - N 
FOR I - 1 TO N 
X - PR(I) : Y - PV(I) 

FOR J - 2 TO 2 * D + 1 
A ( J) - A ( J) + X A (J - 1) 

NEXT J 

FOR K - 1 TO D + 1 

R(K, D + 2) - t(K) + Y * X * (K - 1) 
t(K) - t(K) + Y * X * (K - 1) 

NEXT K 

t (D + 2) - t (D + 2) + Y * 2 
NEXT I 

FOR J - 1 TO D + 1 
FOR K - 1 TO D + 1 
R(J, K) - A ( J + K - 1) 

NEXT K 
NEXT J 

FOR J - 1 TO D + 1 
K - J 

280 IF R (K, J) <> 0 THEN 320 
K - K + 1 

IF K <- D + 1 THEN 280 
PRINT "NO UNIQUE’ SOLUTION" 

GOTO 790 

320 FOR I - 1 TO D + 2 
S - R ( J, I) 

R(J, I) - R(K, I) 

R (K, I) - S 
NEXT I 

Z - 1 / R(J, J) 

FOR I - 1 TO D + 2 
R(J, I) - Z * R ( J, I) 

NEXT I 

FOR K - 1 TO D + 1 
IF K - J THEN 470 
Z - -R(K, J) 

FOR I - 1 TO D + 2 

R(K, I) - R(K, I) + Z * R ( J, I) 

NEXT I 
470 NEXT K 
NEXT J 
PRINT 

PRINT " CONSTANT R(l, D + 2) : CONSTA - R(l, D + 2) : VR(0) 

ONSTA 

FOR J - 1 TO D 

PRINT J; "DEGREE COEFFICIENT R ( J + 1, D + 2) : VR(J) - R(J + 1, D + 2) 

NEXT J 
PRINT 
P ■ 0 

FOR J « 2 TO D + 1 

P - P + R ( J, D + 2) * (t(J) - A ( J) * t(l) / N) 

NEXT J 

q - t (D + 2) - t (1) A 2 / N 
Z - q - P 
I - N - D - 1 
PRINT 
J - P / q 

PRINT "COEFFICIENT OF DETERMINATION (R-'2) - J 
PRINT "COEFFICIENT OF CORRLELATION SQR(ABS(J)) 

PRINT "STANDARD ERROR OF ESTIMATE SQR(ABS(Z / I)) 

790 PRINT "POLYFIT COMPLETED" 

END SUB 
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FORTRAN CODE 


101 



597 

598 

599 

600 
601 
602 

603 

604 

605 

606 

607 

608 

609 

610 
611 
612 

613 

614 

615 

616 

617 

618 

619 

620 
621 
622 

623 

624 

625 

626 

627 

628 

629 

630 

631 

632 

633 

634 

635 

636 

637 

638 

639 

640 

641 

642 

643 


SUBROUTINE POLYFIT (PR, PV, DEGREE, N, VR,VR0) 
IMPLICIT REAL* 16 (A-H,0-Z) 

IMPLICIT INTEGER (I-N) 

INTEGER D , DEGREE, AU 

DIMENSION A(21),R(13, 14), t(14),VR(10),PR(2100),PV(2100) 
D = DEGREE 
AU = 5 

A(l) = QFLOAT( N ) 

DO 10 1 = 1 , N 
X = PR(I) 

Y = PV(I) 

DO 20 J = 2, 2 * D + 1 
A(J) = A(J) + X ** (J - 1) 

20 CONTINUE 
DO 30 K =1,D+1 

R(K, D + 2) = t(K) + Y * X ** (K - 1) 
t(K) = t(K) + Y * X ** (K - 1) 

30 CONTINUE 

t(D +2)=t(D + 2) + Y ** 2. 

10 CONTINUE 
DO 40 J = 1 , D + 1 
DO 50 K = 1 , D + 1 
R(J, K) = A(J + K -1) 

50 CONTINUE 
40 CONTINUE 
DO 100 J = 1, D + 1 
K = J 

280 IF (R(K, J).NE. 0.0) GOTO 320 
K = K+ 1 

IF (K.LE. D + 1 ) GOTO 280 
WRITE(AU,7) 

7 FORMAT (’ NO UNIQUE SOLUTION’) 

GOTO 790 

320 DO 60 1 = 1 , D + 2 
S = R(J, I) 

R(J,I) = R(K, I) 

R(K, I) = S 
60 CONTINUE 
Z= 1.0 /R(J, J) 

DO 70 1 = 1 , D + 2 
R(JJ) = Z*R(J,I) 

70 CONTINUE 

DO 470 K = 1 , D + 1 
IF (K.EQ.J) GOTO 470 
Z = -R(KJ) 

DO 90 I = 1, D + 2 

R(K, I) = R(K, I) + Z * R(J, I) 
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644 90 CONTINUE 

645 470 CONTINUE 

646 100 CONTINUE 

647 AU = 1 

648 WRITE(AU,17) R(l, D+2) 

649 17 FORMAT (’0 CONSTANT =\E15.8) 

650 CONSTA = R(l, D + 2) 

651 VR0 = CONSTA 

652 DO 1 10 J = 1 , D 

653 WRITE(AU,27) J, R(J+1 ,D+2) 

654 27 FORMAT(lXJ4, ’DEGREE COEFFICIENT =’,E15.8) 

655 VR(J) = R( J + 1 , D + 2) 

656 110 CONTINUE 

657 P = 0 

658 DO 1 20 J = 2 ,D + 1 

659 P = P + R(J, D + 2) * (t(J) - A(J) * t(l) / QFLOAT( N ) ) 

660 120 CONTINUE 

661 q = t(D + 2) - t(l) ** 2. / QFLOAT( N ) 

662 Z = q - P 

663 I = N - D - 1 

664 PJ=P/q 

665 WRJTE(AU,37) PJ 

666 37 FORMAT (’ COEFFICIENT OF DETERMINATION (R*2)=\E15.8 ) 

667 SQAJ = QSQRT( QABS( PJ ) ) 

668 SQAZI = QSQRT(QABS(Z/QFLOAT(I))) 

669 WRJTE(AU,47) SQAJ 

670 47 FORMAT (’ COEFFICIENT OF CORRELATION =\E15.8) 

671 WRITE(AU,57) SQAZI 

672 57 FORMAT (’ STANDARD ERROR OF ESTIMATE =\E15.8) 

673 790 PRINT VPOLYFIT COMPLETED’ 

674 RETURN 

675 END 
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RUK and RUK4 Subroutine Descriptions 


The RUK and RUK4 subroutine are fourth order Runge-Kutta integration routines. They call D 
ER1 and DER, respectively, and integrate the equations of motion of the spacecraft. 


RUK and RUK4 Subroutine Passed Variables 

SUB RUK (X(), N, dt) STATIC 
SUB RUK4 (X(), N, dt) STATIC 


dt Integration step size (seconds). 

NUM Order of the X state vector for the Runge-Kutta. 

X() State Vector of the spacecraft’s position and velocity (km and km/s). 
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BASIC CODE 
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DEFDBL A-H, J-M, O-Z 

SUB R0K4 (X(), N, dt) STATIC 

DIM D (6) , F (6) , 0(6), DX(6) 

CALL DER(X() , D() ) 

FOR I - 1 TO N 

D (I) - D (I) * dt: 0(1) - X(I) + .5# * D(I) 

NEXT I 

CALL DER (O () , F() ) 

FOR I - 1 TO N 

F (I) - F (I) * dt: D(I) - D(I) + 2# * F(I): 0(1) - X(I) + .5# * F (I) 

NEXT I 

CALL DER (O (} , F() ) 

FOR I - 1 TO N 

F (I) - F (I) * dt: D (I) - D (I) + 2# * F(I) : 0(1) - X(I) + F(I) 

NEXT I 

CALL DER (O ( ) , F() ) 

FOR I - 1 TO N 

X(I) - X(I) + (D (I) + F (I) * dt) / 6# 

NEXT I 
END SOB 
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706 SUBROUTINE RUK4 (X, N, dt, MU, ACCEL 1 , DM, DE,MUM, 

707 * MUE,TH,VR,thrst,CJ, ACCEL, RANGE, DIRECT, ALPHA, 

708 * RANGEOFF,theta,RJAC,VTN,DIST,VRR,VRN,VRO,TT,VT) 

709 IMPLICIT REAL* 1 6 (A-Z) 

710 INTEGER I,N,TH,thrst,VTN,VRN 

711 DIMENSION D(4),F(4),U(4),X(4),VR(10),VT(5),RJAC(3) 

712 CALL DER(X, D,DM,DEMUM,MUE,TH,VR,thrst,CJ,ACCELJLANGE, 

713 * DIRECT ,ALPHA,RANGEOFF, theta, RJ AC, VTN,VT,DIST,VRR,VRN,VRO,TT) 

714 DO 10 1 = 1 , N 

715 D(I) = D(I) * dt 

716 U(I) = X(I) + .5 * D( I ) 

717 10 CONTINUE 

7 1 8 CALL DER(U, F,DM,DE,MUM,MUE,TH,VR,thrst,CJ,ACCEL,RANGE, 

719 * DIRECT, ALPHA, RANGEOFF,theta,RJ AC, VTN,VT,DIST,VRR,VRN,VRO,TT) 

720 DO 20 I = 1 , N 

721 F(I) = F(I) * dt 

722 D(I) = D(I) + 2. * F(I) 

723 U(I) = X(I) + .5 * F(I) 

724 20 CONTINUE 

725 CALL DER(U, F,DM,DE,MUM,MUE,TH,VR,thrst,CJ,ACCEL,RANGE, 

726 * DIRECT, ALPHA,RANGEOFF, theta, RJ AC, VTN,VT,DIST,VRR,VRN,VR0,TT) 

727 DO 30 1 = 1 , N 

728 F(I) = F(I) * dt 

729 D(I) = D(I) + 2. * F(I) 

730 U(I) = X(I) + F(I) 

731 30 CONTINUE 

732 CALL DER(U, F,DM,DE,MUM,MUE,TH,VR,thrst,CJ,ACCELJLANGE, 

733 * DIRECT, ALPHA,RANGEOFF,theta,RJAC,VTN,VT,DIST,VRR,VRN,VR0,TT) 

734 DO 401= 1.N 

735 Xa> =X(I) + (D(I) + F(I) *dt)/6. 

736 40 CONTINUE 

737 RETURN 

738 END 
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SPIRAL Subroutine Description 


The subroutine SPIRAL is used to generate the parametric curves of the radial and tangential 
velocity components for the reference capture spiral. The direction of the cislunar trajectory 
determines the governing gravitational parameter, MU, for the spiral trajectory generation. The 
estimated final mass about the target planet is taken to be 80 % of the spacecraft’s chosen initial 
mass. The subroutine begins an integration loop using a negative mass flow for the spacecraft’s 
propulsion system. A spiral trajectory is generated out from the target planet of the spacecraft 
gaining mass as it goes. This mimics the ideal spiral capture with the spacecraft losing mass as it 
settles into the capture orbit. The integration routine is a fourth order Runge-Kutta that calls the 
derivative subroutine, DERI. DERI contains the equations of motion for a two-body trajectory. 
No perturbational forces, other than the spacecraft’s acceleration, are included in the two- 
dimensional equations of the motion. SPIRAL calculates the acceleration level, the two-body 
energy, and the tangential and radial velocity components. These velocity components are 
determined by finding the angle, theta, between the radial vector and the velocity vector using 
the dot products rule, 

r x v = 1 r | |v| cos (theta) 

where r is the radius, and v is the velocity. The radial component of the velocity is found by 
multiplying the cosine of theta by VMAG, the magnitude of the velocity vector. The tangential 
component of the velocity vector is similarly found by multiplying VMAG by the sin of theta. 
Every tenth integration the velocity components and radial distance is captured into three arrays, 
PVR(), PVT(), and PR(). When the two-body energy of the spiral trajectory is non-negative the 
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spacecraft’s thrust is turned off and the spacecraft is assumed to be following a parabolic path. 
The integration continues until the spacecraft passes the range flag, RANGESC. POLYFIT 
develops a polynomial curve fit of the radial velocity component. The tangential velocity 
component is fit to either a power, exponential, logarithmic, or linear curve in the subroutine 
CURVE depending on which equation best fits the data. When the parametric curve fitting is 
complete the program records the velocity curve coefficients into VR() and VT() and returns 
control to IO. 
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SPIRAL Subroutine Internal and Shared Variables 


SUB SPIRAL (X(), DIRECTION, RANGESC, VT(), VR() 

ctheta Cosine of the angle between the radius and velocity vector. 

DIRECTION Numeric indicator of the direction of the trajectory generation, Earth to Moon (0) 
or Moon to Earth (1). 

dt Integration step size (seconds). 

e Two-body energy, sum of the potential and kinetic energy. 

J The array size counter for the radial and tangential velocity, and the radial 

distance arrays. 

PR() The array containing the radial distance from the capture planet (km/1000). 

PVR() The array containing the radial component of velocity, vrad, during the reverse 

integration spiral (km/s). 

PVT() The array containing the tangential component of velocity, vtan, during the 

reverse integration spiral (km/s). 

RANGESC Range from the central body for the reverse integration spiral to be generated 
before stopping (km). 

RDOTV Radial vector dot multiplied with the velocity vector. 

RMAG Magnitude of the spacecraft’s radius vector from the capture planet (km). 

SCMASS1 Starting mass for reverse integration of spacecraft from target planet, 80% of 
chosen initial mass. 

TT Total time of trajectory generation (seconds). 

VMAG Magnitude of the spacecraft’s velocity vector (km/s). 

VMAG2 Square of the velocity magnitude of the spacecraft in relation to the capture planet 
(km/s) A 2. 

VR() The array of the parameterized radial velocity coefficients. 

vrad Component of the spacecraft’s velocity in the radial direction (km/s). 

VT() The array of parameterized tangential velocity coefficients. 
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vtan 

vthet 

w 

X() 


Component of the spacecraft’s velocity in the tangential direction (km/s). 
Angle between the radius and velocity vector (radians). 

Counter for sampling the capture spiral trajectory. 

State Vector of the spacecraft’s position and velocity in the rotating x,y 
nates (km and km/s). 


coordi- 
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SPIRAL Subroutine Flowchart 


SPIRAL 



Set gravitational parametei 
MU, for reference spiral 



Call POLYFIT 


Call CURVE 


RETURN 
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BASIC CODE 
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DEFDBL A-Z 

SUB SPIRAL (X(), DIRECTION, RANGES C, VT(), VR()) 

REDIM PR (2100) , PVR(2100), PVT(2100) 

IF DIRECTION - 1 THEN 
MU - MUE 
ELSE 

MU - MUM 

END IF 

I . 0: w - 5: J-0: dt-3 

SCMASS1 - SCMASS * .8# 'total -s/c mass (kg) 

SCREEN 2: WINDOW (70000, 50000 )- (-70000, -50000) 

CIRCLE (0, 0) , X (2) , • 1 
DO 

RUK X(), NUM, dt 

RMAG - SQR (X (1) A 2 + X(2) A 2) 

dt - .02# * RMAG 

TT - TT + dt: I - I + 1 

VMAG2 - X (3) A 2 + X(4) A 2 

e - VMAG2 / 2# - MU / RMAG 

ACCEL1 - THRUST / (SCMASS1 + Mdot * TT) 'acceleration of the s/c (km/s A 2) 
IF (e >- 0#) THEN ACCEL1 - 0# 

IF w - 10 THEN 

VMAG - SQR (VMAG2) 

RDOTV - X(l) * X (3) + X (2) * X(4) 

Ctheta - RDOTV / (RMAG * VMAG) 
vthet - ACOS (ctheta) 
vrad - VMAG * COS (vthet) 
vtan - VMAG * SIN (vthet) 

J - J + 1: SHOW X() 

PR ( J) - RMAG / 1000#: PVR(J) - vrad: PVT(J) - vtan: w - 0 
END IF 
w ■■ w + 1 

LOOP WHILE RMAG < RANGESC 

PRINT " done”, RMAG, vrad, vtan, I, J, e 

POLYFIT PRO, PVR ( ) , VRN, J, VR ( ) 

PRINT "The coefficients of vrad vs rad polynomial” 

CURVE PRO, PVT ( ) , J, VT() 

END SUB 
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FORTRAN CODE 
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740 SUBROUTINE SPIRAL (X, DIRECT, RANGESC, VT, VR, MUE, 

-741 + MUM, SCMASS, dt, THRUST, TT, Mdot, 

742 + MU, ACCEL1,NUM,DM,DE,VRN,VTN,VR0) 

743 IMPLICIT REAL* 1 6 ( A-H,0-Z) 

744 REAL* 1 6 MU,MUM,MUE,Mdot 

745 INTEGER W,VRN, DIRECT, AU.VTN 

746 DIMENSION X(4), PR(2100), PVR(2100), PVT(2100), VR(10),VT(5) 

-747 DIMENSION ABCD(IOO) 

748 AU = 5 

749 IF (DIRECT.EQ. 1 ) THEN 

—750 MU = MUE 

751 ELSE 

752 MU = MUM 

-753 ENDIF 

754 w = 5 

755 J = 0 

_756 dt = 3.0 

757 SCMASS 1 = SCMASS * .8 

758 DO WHILE ( RMAG .LT. RANGESC ) 

759 CALL RUK (X,NUM, dt,MU,ACCELl) 

760 RMAG = QSQRT(X( 1 ) ** 2. + X(2) ** 2.) 

761 dt = .02 * RMAG 

762 TT = TT + dt 

763 VMAG2 = X(3) ** 2. + X(4) ** 2. 

764 e = VMAG2 / 2. - MU / RMAG 

765 ACCEL 1 = THRUST/(SCMASS 1 + Mdot * TT) 

766 IF (e.GE.0.0) ACCEL1 =0.0 

767 IF (w.EQ.10) THEN 

768 VMAG = QSQRT(VMAG2) 

769 RDOTV = X( 1 ) * X(3) + X(2) * X(4) 

770 ctheta = RDOTV / (RMAG * VMAG) 

77 1 vthet = QACOS(ctheta) 

~772 vrad = VMAG * QCOS(vthet) 

773 vtan = VMAG * QSIN(vthet) 

774 J=J+1 

775 PR(J)=RMAG/1 000.0 

776 PVR(J) = vrad 

777 PVT(J) = vtan 

-778 w=0 

779 ENDIF 

780 w= w+ 1 

-781 END DO 

782 CALL POLYFIT (PR,PVR,VRN,J,VR,VR0) 

783 AU = 1 

-784 WRTTE(AU,107) 

785 107 FORMAT (TThe coefficients of vrad vs rad polynomial’ ) 

786 CALL CURVE(PR,PVT J,VT,VR0,VTN) 
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787 RETURN 

788 END 



